[gmx-users] Bonds of chloroform molecules messed up after EM

[Justin Lemkul]

On 10/4/18 6:38 PM, Lili Zhang wrote:
> Hello Gromacs users,
>
> I am trying to run MD  simulation of chloroform in hydrated POPC in
> CHARMM36 ff. I did not observe anything wrong with the bond structure of
> chloroform after insert molecules or solvation. But after EM, when I looked
> em.gro in VMD, there were Cl-Cl bonds or Cl-H bonds formed on some
> chloroform molecules, which had not been there for npt.gro. I assume I need
> to use SHAKE to constrain those bonds during EM, but no completely sure.
>
> And I also want check if my "Adding a residue to a force field" process is
> done correctly. The following is what I did.
>
> ========================================================================
> The chloroform atomtypes CCM HCM and CLCM are not included in the ff, so I
> first added them to the atomtypes.atp as
>    CLCM     35.45000 ;
>     HCM      1.00800 ;
>     CCM     12.01100 ;
>
> Second, I added LJ parameters to ffnonbonded.itp as
> ;type atnum         mass   charge ptype           sigma  epsilon
> CLCM     17    35.450000    0.000  A  0.172067736185  1.25520
>   CCM      6    12.011000    0.000  A  0.170071301248  0.42677
>   HCM      1     1.008000    0.000  A  0.110044563280  0.08368
>
> Third, I added bond info to ffbonded.itp as
> [ bondtypes ]
> ;      i        j  func           b0           kb
>       CCM     CLCM     1   0.17580000    418400.00
>      CCM       HCM     1   0.11000000    418400.00
>
> [ angletypes ]
> ;      i        j        k  func       theta0       ktheta          ub0
>       kub
>      CLCM      CCM     CLCM     5   111.300000     0.000000   0.00000000
>      0.00
>       HCM      CCM     CLCM     5   107.570000     0.000000   0.00000000
>      0.00

The force constants on your angles are zero, so you effectively have no 
angle potential and your molecule will surely distort as a consequence.

If this is the Dietz and Heininger model for CHCl3, it relies on a rigid 
geometry that includes constrained bonds between all atoms. As I recall 
we have discussed this before, it may not work in GROMACS. But if you're 
going to convert the topology to a "normal" approach of just covalent 
bonds and valence angles, you have to actually parametrize the bonds and 
angles, not just assign arbitrary values (which looks like the case for 
kb in the bonds, as well).

-Justin

-- 
==================================================

Justin A. Lemkul, Ph.D.
Assistant Professor
Virginia Tech Department of Biochemistry

303 Engel Hall
340 West Campus Dr.
Blacksburg, VA 24061

jalemkul at vt.edu | (540) 231-3129
http://www.thelemkullab.com

==================================================