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 BOXTYPE = RECT
SET XTLTYPE = TETRAGONAL
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 NLIPTOP = 19
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)
SET NWATER = 924
SET ISGLIP = NO
SET ISMEMB = YES


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.
******************************************
BOMLEV ( 0) IS REACHED - TERMINATING. WRNLEV IS 5


*** 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.

Last edited by ASmittie; 03/14/19 08:08 AM.