Previous Thread
Next Thread
Print Thread
Page 1 of 2 1 2
Joined: Mar 2019
Posts: 14
A
Forum Member
OP Offline
Forum Member
A
Joined: Mar 2019
Posts: 14
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.
Joined: Sep 2003
Posts: 8,605
Likes: 23
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,605
Likes: 23
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.

Last edited by rmv; 03/15/19 04:20 PM. Reason: chirality

Rick Venable
computational chemist

Joined: Mar 2019
Posts: 14
A
Forum Member
OP Offline
Forum Member
A
Joined: Mar 2019
Posts: 14
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:

SELECTED IMAGES ATOMS BEING CENTERED ABOUT 0.000000 0.000000 0.000000
UPDECI: Nonbond update at step 2553
Using Image CUBE search
TOTAL ENERGY CHANGE EXCEEDED

20. KCAL AND 10% OF THE TOTAL KINETIC ENERGY IN THE LAST STEP
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).
     NOTE!! THIS FILE C A N N O T BE USED TO RESTART A RUN!!!



   ***** LEVEL -2 WARNING FROM *****
   ***** ENERGY CHANGE TOLERANCE EXCEEDED
   ******************************************
   BOMLEV ( 0) IS REACHED - TERMINATING. WRNLEV IS 5




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

INTERNAL COORDINATES
TITLE> * GENERATED BY CHARMM-GUI V2.0 ON DEC, 19. 2018. JOB
TITLE> * INPUT FILES FOR EQUILIBRATION
TITLE> *

NUMBER OF INTERNAL COORDINATES : 0

  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
SELECTED COORDINATES COPIED TO THE COMPARISON SET.

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

Joined: Sep 2003
Posts: 8,605
Likes: 23
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,605
Likes: 23
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


Rick Venable
computational chemist

Joined: Mar 2019
Posts: 14
A
Forum Member
OP Offline
Forum Member
A
Joined: Mar 2019
Posts: 14
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.

Joined: Sep 2003
Posts: 8,605
Likes: 23
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,605
Likes: 23
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.


Rick Venable
computational chemist

Joined: Mar 2019
Posts: 14
A
Forum Member
OP Offline
Forum Member
A
Joined: Mar 2019
Posts: 14
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?

Joined: Sep 2003
Posts: 8,605
Likes: 23
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,605
Likes: 23
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 charmm.org 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


Rick Venable
computational chemist

Joined: Sep 2003
Posts: 4,850
Likes: 6
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 4,850
Likes: 6
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.


Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden
Joined: Sep 2003
Posts: 8,605
Likes: 23
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,605
Likes: 23
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.


Rick Venable
computational chemist

Page 1 of 2 1 2

Moderated by  BRBrooks, lennart, rmv 

Link Copied to Clipboard
Powered by UBB.threads™ PHP Forum Software 7.7.5
(Release build 20201027)
Responsive Width:

PHP: 7.3.31-1~deb10u1 Page Time: 0.008s Queries: 35 (0.003s) Memory: 0.7948 MB (Peak: 0.9003 MB) Data Comp: Off Server Time: 2022-09-27 01:36:05 UTC
Valid HTML 5 and Valid CSS