[gmx-users] Acetonitrile with CHARMM ff

[Sonia Milena Aguilera Segura]

Dear Justin, 

Thank you very much for your answer. I tried both approaches: erasing the dihedral from the topology (first approach) and modifying the ffbonded.itp as described in the previous mail, without adding any improper (second approach). I ran a minimization, 200 ps NVT, 200 ps NPT, and 1 ns MD. In general the results seem pretty much similar with respect to the energy terms:
   
                           1st             2nd
E.Pot			-22425.4	-22452.9
E. Kin			6940.56	         6928.76
E. Tot			-15484.8	-15524.2
Temp (K)	        301.302	           300.79
Pressure (bar)	        2.35245	        -1.23941
Density			788.145 	788.793
diffusivity x10-5 cm2/s	2.0584 (+/- 0.0720)	2.1406 (+/- 0.1427)
dielec. Constant	20.6455	        19.8552

However, I compared the density and it's higher with respecter to the experiment. Also, I noticed that the values T are slightly higher and P seems to show really different behavior. I checked the results for the NPT equilibration and the values were 8.8 and 14.5, respectively. The temperatures during the NVT were  297.565	and 297.799. Moreover, dielectric constants are underestimated. Can it be that I am using the wrong parameters?  Here is my mdp file for the MD (same for NVT and NPT, adjusting the coupling terms):

; Run control
integrator               = sd       ; Langevin dynamics
tinit                    = 0
dt                       = 0.0005
nsteps                   = 40000000   ; 20 ns
nstcomm                  = 100
; Neighborsearching and short-range nonbonded interactions
cutoff-scheme            = verlet
nstlist                  = 20
ns_type                  = grid
pbc                      = xyz
rlist                    = 1.2
; Electrostatics
coulombtype              = PME
rcoulomb                 = 1.2
; van der Waals
vdwtype                  = cutoff
vdw-modifier             = potential-switch
rvdw-switch              = 1.0
rvdw                     = 1.2
; Apply long range dispersion corrections for Energy and Pressure
DispCorr                  = EnerPres
; Spacing for the PME/PPPM FFT grid
fourierspacing           = 0.12
; EWALD/PME/PPPM parameters
pme_order                = 6
ewald_rtol               = 1e-06
epsilon_surface          = 0
; Temperature coupling
; tcoupl is implicitly handled by the sd integrator
tc_grps                  = system
tau_t                    = 1.0
ref_t                    = 298.15
; Pressure coupling is on for NPT
Pcoupl                   = Parrinello-Rahman 
tau_p                    = 1.0
compressibility          = 4.5e-05
ref_p                    = 1.0 
; Do not generate velocities
gen_vel                  = no 
; options for bonds
constraints              = none  ; we only have C-H bonds here
; Type of constraint algorithm
constraint-algorithm     = lincs
; Constrain the starting configuration
; since we are continuing from NPT
continuation             = yes 
; Highest order in the expansion of the constraint coupling matrix
lincs-order              = 12

Than you very much for any advice, 

Best regards, 

Sonia Aguilera
PhD Student
ENSCM
> 
> On 7/6/17 6:34 AM, Sonia Milena Aguilera Segura wrote:
> > Dear Justin,
> > 
> > Thank you for your reply. I went trough the CHARMM forum and I found this
> > post of 2007,
> > https://www.charmm.org/ubbthreads/ubbthreads.php?ubb=showflat&Number=28545,
> > which basically underlines the same issues I am having. I read this post
> > before posting here, and this why I came with the idea of modifying the
> > dihedral ff parameters, as suggested in the forum. Apparently, it may work
> > if I set the amplitude k=0 and I distort the linear angle a little (maybe
> > 179.9) in the initial structure. The thing is, as I say before, I have no
> > experience modifying force fields, so I wouldn't what to put and where. My
> > guess is that I would need to go the ffbonded.itp file and modify and add
> > the following lines in this way:
> > 
> > [ angletypes ]
> > ;      i        j        k  func       theta0       ktheta          ub0
> > kub
> > CG331    CG1N1    NG1T1     5   179.900000   177.401600   0.00000000
> > 0.00
> > 
> > [ dihedraltypes ]
> > ;      i        j        k        l  func         phi0         kphi  mult
> >       HGA3     CG331    CG1N1  NG1T1     9     0.000000     0.000000     1
> >    

> 
> You can modify the [angletype] parameters.  A simple energy minimization will
> correct the structure.  With respect to the dihedrals, as noted in that link,
> they simply shouldn't be there.  You can do as above with a k=0 parameter
> (don't
> add an improper, though harmless in this case it is incorrect to add these
> terms
> to methyl groups).  Or, if you want to avoid adding dummy parameters, just
> remove the dihedral lines from your .top file altogether.
> 
> -Justin
>