CHARMM Development Project
I am a grad student in a quantum computational research lab at University. Our research has led me into the realm of classical dynamics in order to support findings from gas-phase calculations. We are hoping to determine explicit solvent effects on conformational dynamics of a small, biologically active organic molecule consisting of a ring structure and a short hydrocarbon tail. I have built a hydrated DPPC membrane using the CHARMM-gui bilayer builder.

This is my crystal structure as found in step5_assembly.str:
SET A = 34.6285541
SET B = 34.6285541
SET C = 69
SET ALPHA = 90.0
SET BETA = 90.0
SET GAMMA = 90.0
SET ZCEN = 0.0
SET NLIPBOT = 18 <<< this leaf has my target organic molecule centered in the xy-plane (z-position not set at equilibrium distance from bilayer center)

I want to keep the conformation of my target molecule fixed until the equilibration steps are complete and it has diffused to an equilibrium position in the bi-layer. I am looking to apply an absolute constraint to the internal coordinates for my target molecule while allowing for whole-molecule translation/rotation.

So far I have tried the following approaches: (I have added portions of each step6.1_equilibration.inp file relating to the segment definitions and constraints.)

*** constrain conformationally significant dihedral angles with CONS DIHE and a large force constant.

define DIHA sele ( segid HETA .and. ( type C23 .or. type C19 .or. type C17 .or. type C15 ) ) end
define DIHB sele ( segid HETA .and. ( type C12 .or. type C9  .or. type C15 .or. type C17 ) ) end

cons dihe force 1000.0 sele DIHA end
cons dihe force 1000.0 sele DIHB end

--> this method still resulted in deviation of specified dihedrals from ideal angle values. What sort of values for force constant are too large? I'm very cautious to apply large force constants during equilibration dynamics.

*** load previously defined coords into COMParison set; build IC table from COMP; apply ic cons with high force constants for BOND, ANGL, DIHE, IMPR, and a large EXPO value.

define RING sele segid HETA .and. -
( ( type C1   .or. type C2   .or. type C3   .or. type C4   .or. -
    type C5   .or. type C6   .or. type C8   .or. type C9   .or. -
    type C10 .or. type C12 .or. type C15 .or. type C17 .or. -
    type C19 .or. type C20 .or. type C23 .or. type C24 .or. -
    type C25 .or. type C26 .or. type C27 ) .or. ( type O ) .or. -
  ( type H1   .or. type H2   .or. type H3   .or. type H4   .or. -
    type H5   .or. type H6   .or. type H7   .or. type H8   .or. -
    type H10 .or. type H11 .or. type H12 .or. type H13 .or. -
    type H14 .or. type H17 .or. type H23 .or. type H26 .or. -
    type H28 .or. type H29 .or. type H36 .or. type H37 .or. -
    type H38 .or. type H39 .or. type H40 .or. type H41 .or. -
    type H42 .or. type H43 .or. type H44 ) ) end

coor copy sele RING comp end
ic fill comp
cons ic bond 1000 angl 1000 dihe 1000 impr 1000 expo 10

--> I keep encountering the null set error when I attempt to build the IC table from comp. Is there a direct method that I can build the IC table from only a defined set of atoms?

I tried the solution offered by rmv here but the command:
gener a setup
CHARMM> gener a setup
GENPSF> Segment 3 has been generated. Its identifier is TIP3.

***** LEVEL 0 WARNING FROM *****
***** PSF modified with image centering active: Now off.

*** isolate the target molecule with CONS DROP, although I suspect that this is not the correct use of that constraint.

Please let me know what additional information (if any) is needed to clarify my goal.
The bilayer system seems a bit too small, and may lead to periodicity artifacts.

It's not clear if you have an IC table established. How did you get one for your small molecule?

I would not apply restraints to bonds and angles with CONS IC, and suggest starting with the default exponent. You may not need IMPR either, unless the molecule has a carbonyl group.

Use CRYSTAL FREE before GENERATE, then re-establish the CRYSTAL setup for PBC afterwards.

P.S.: It occurred to me that maintaining a chiral center might be a good reason to restrain selected IMPR IC entries.
I understand that my system is small and that this may produce artifacts. Ideally if I can get the desired effects from my constraints I will embed the molecule in a much larger system. Unless you think that starting out in a larger box will avoid potential force/energy issues. Such as in my original post I described one approach by defining two dihedral angles, DIHA & DIHB, with a force of 1000.0. And even a force of 4000.0 kJ/mol gives deviation of approx. +-5 degrees for each dihedral. I found that if I with a force of 5000.0 the calculation crashes within the first equilibration step with this error:

UPDECI: Nonbond update at step 2553
Using Image CUBE search

PREVIOUS E = -914.4 CURRENT E = -174.1 KINETIC = 6861.
WRIDYN: RESTart file was written at step 2573

     Writing RESTART FILE with previous and current coordinates,
     which may be read by: READ COOR DYNR .... (see io.doc).

   ***** LEVEL -2 WARNING FROM *****

As for the IC table:

I do not have a pre-constructed IC table for my target molecule. Also, my IC table is empty from my initial inputs in step6.1_equilibration.inp.

CHARMM> print ic



  N   I   J   K   L   R(I(J/K))   T(I(JK/KJ))   PHI   T(JKL)   R(KL)

I tried using CRYSTAL FREE followed by IC GENERate COMP and was able to successfully fill the IC table from the comparison set.

I noticed, however, that my command to save into comparison set only those atoms that comprise my target molecule was not performing as I had intended. With the assumption that the comp set would be overwritten, I was using:

coor copy sele segid HETA end comp
CHARMM> coor copy sele segid HETA end comp
SELRPN> 72 atoms have been selected out of 7654

This is the last line from PRINT COOR COMP after issuing the above command.
7654 962 TIP3 H2 9999.0000000000 9999.0000000000 9999.0000000000 TIP3 924 0.0000000000

So clearly I am not replacing the full comparison set with only a selection of atoms. Am I missing something? Can the comparison set be overwritten as I want?

Is this even a reasonable route for accomplishing my goal? (Which is to prevent any conformational changes to my target molecule during equilibration steps while still allowing it to translate/rotate within the bilayer.)
Changes of 10 or 20 deg for a dihedral are not conformational changes in most cases, but more in the line of expected fluctuations. As you've observed, overly aggressive restraints can lead to problems. BTW, the energy units in CHARMM are kcal/mol, not kJ/mol, and the force constant is per radian, and not per degree. I would not restrain the hydrocarbon tail at all.

Both MAIN and COMP coordinate sets both contain all atoms defined in the PSF; the atom selection with COOR COPY limits what is copied, but not which atoms are present in the result.

Note that the "intcor" doc file warns that IC GENERATE can have issue with ring systems. IMHO, the best IC tables are hand edited, such as those found in the protein and nucleic acid topology files.

The system size is more of an issue for any attempt to estimate lateral diffusion in a bilayer:

"Lipid and Peptide Diffusion in Bilayers: The Saffman-Delbrück Model and Periodic Boundary Conditions," R. M. Venable, H. I. Ingólfsson, M. G. Lerner, B. S. Perrin, B. A. Camley, S. J. Marrink, F. L. H. Brown and R. W. Pastor, The Journal of Physical Chemistry B, 121 (15) pp. 3443-3457 (2017). DOI
More precisely put, I want to prevent any changes to bond distances and angles for the target molecule until the equilibration steps have finished and I am into production steps.

I mention diffusion because I don't currently know its equilibrium position along the bilayer normal. This is relevant because this molecule has a polar functional group opposite from a hydrocarbon chain.
Placing the molecule with the hydrocarbon chain embedded in the lipid chains and the polar part in the lipid head group region would be a reasonable choice. Umbrella sampling is often used to determine a free energy z profile to refine the optimal position.

For a well built model, restraints on bonds and angles should not be needed; they are fairly rigid, and will quickly adopt the FF equilibrium values in the absence of any restraints.

A concern for model building with ring systems is checking for penetration of rings by atoms from another molecule placed in the same vicinity.
The atomic coordinates for the target molecule were taken from a conical intersection of a gas-phase quantum calculation. We want to embed the molecule in a membrane and study the relaxation pathway and compare to the gas-phase dynamics.

I am looking for a means to prevent any changes to the internal coordinates of the molecule during equilibration. Is the only method for doing this by hand-building the IC table? And can I place an absolute constraint on that IC table?
Most of the "CONS" commands are actually restraints, and not formal constraints; they add a force, rather than remove degrees of freedom.

In principle one could use the SHAKE constraint, but I believe that it is not respected for heavy atom bonds in some cases. I'm fairly certain it does not work with the DOMDEC engine, and it may even be in an issue in general for parallel energy evaluations. (I haven't tried or investigated that for myself.)

There is also the SHAPES facility (which I have also never used), but it appears that doc file has issues that the automated HTML generator used at can't cope with. You can browse the file in the CHARMM distribution at doc/shapes.doc, which is a text file (not a Word file). A version of the page can also be found here
SHAKE BONDs is applied to all bonds, and there is also SHAKE angles, which is not generallly recommended since angles and dihedrals are coupled. It might be OK just for setting things up though.
With the use of the appropriate atom selections, you can limit SHAKE BONDS to only your target molecule.

Note that SHAKE FAST only supports bonds involving H atoms; since DOMDEC requires its use, it appears that SHAKE BONDS will not work with that fast MD engine.
I have been trying to use shake to constrain the bond distances of my target molecule, however this does not prevent the bond angles from changing. When I try shake angle then all the atoms in the molecule move wildly.

It is baffling me that I am unable to find a function in CHARMM to perform a task that is not uncommon (preventing any changes to a molecule's shape).

Am I looking at the wrong commands? I have tried:

cons dihe
cons harm abso
cons IC

shake bond
shake angl
Which biomolecular MD programs have a function to prevent any changes to a molecule's shape, and how do these function differ from the constraints and restraints available in CHARMM?

You have to remember that the main purpose of MD is to change the conformation...

All atoms in the molecule should not move wildly when you SHAKE ANGLES if the angles in the structure are close to what SHAKE tries to impose; I have however never used SHAKE ANGLES myself.
The task is common, but not to the tolerances you seem to require. The concerns are more typically preserving dihedral states and chirality with restraints during equilibration, and relying on the balls and springs force field for the rest.

One option might be defining custom atoms for your target molecule, so that you can increase the force constants for the bond, angle, and dihedral terms to achieve the degree of stiffness you want. However, this would require a non-trivial effort.

As far as I know, only the SHAPES facility was intended for rigid body dynamics, but I don't know that much about it.

If you wish to allow whole molecule translation and rotation, CONS HARM BESTFIT might be an option instead of CONS HARM ABSO, but I doubt it would achieve the tolerances you want.
GROMACS has this function, however that software is designed for proteins. TurboMole has this function, but is not a biomolecular MD suite.

Lennart, that wasn't a very helpful remark. I am aware that the purpose of MD is to change the conformation... there are also times when you want to change only a portion of the conformation and keep another portion rigid.

rmv: I will try SHAPEs. I checked out that link you provided and am still figuring out the syntax for my needs. I have already tried CONS HARM BESTFIT and didn't get the desired effect.
CHARMM was also designed for biopolymers, (years before GROMACS), but rigid body dynamics has not been pursued very much.
Instead of trying to fix the molecule in place during step6_equilibration, could I let those steps run unconstrained, do one step in production and then edit the coordinates/velocities from my gas phase dynamics into the restart file? I realize that this will require some efforts with positioning in the membrane.

My other idea is read the initial coord/velocities into comparison, do the equilibration and then read the comp set back into the system. Would this work?
© CHARMM forums