Re: [gmx-users] Acetonitrile with CHARMM ff

[Justin Lemkul]

On 7/13/17 9:16 AM, Sonia Milena Aguilera Segura wrote:

Message: 2
Date: Thu, 13 Jul 2017 08:56:14 -0400
From: Justin Lemkul 
To: Discussion list for GROMACS users 
Subject: Re: [gmx-users] umbrella sampling
Message-ID: <995aa507-8875-8c41-674b-84eded63e...@vt.edu>
Content-Type: text/plain; charset=utf-8; format=flowed




Dear Justin,

I increased the size of the simulation box to 4 nm. Indeed, the values of
pressure improved (averages around -1 bar or 1.5, 2 bar, or so). However,
the T keeps being overestimated (302 K). During NVT I got the right value,
so I decided to run the 200 ns NPT equilibration with Berendsen barostat
instead of P-R. I got the right T value, around 298.3 and a pressure of
1.5 more or less. Then I launched a continuation test with 200 ps MD with
P-R (I couldn't use the -e option because it gives me the error Could not
find energy term named 'Box-Vel-XX', moreover, I didn't use P-R before so
I guess I shouldn't expect stored values for it?. But I am using -t ).
Despite I got low pressure averages around 0.5 bar, the temperature raised
again to 301-302. This happens very early in the simulation, which seems
to indicate that for sd integrator/thermostat and P-R barostat there is
something going on. I am getting practically the same density in between
simulation, so I guess I can say t

  ha

   t in term of P, the system has been equilibrated. However, what can I do
   to get the right T for this system? If I was able to get the right T
   with Berendsen barostat, I don't understand what's wrong when I change
   to P-R.



Sounds buggy.  What version of GROMACS are you using?  There were temperature
issues with the Langevin integrator in previous versions, but those should
have
been long since solved.


I am using version 5.1.2. I checked the release notes for version 5.1.3 and 5.1.4 
and the only thing related with sd integrator (none) or P-R barostat is this one 
'Fixed Parrinello-Rahman with nstpcouple > 1'. Can this be the cause?



I don't know.  Upgrade to 2016.3 and try again.

-Justin

--
==

Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 629
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalem...@outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul

==
--
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Acetonitrile with CHARMM ff

[Sonia Milena Aguilera Segura]

- Original Message -
> From: "gromacs org gmx-users-request" 
> 
> To: "gromacs org gmx-users" 
> Sent: Thursday, July 13, 2017 2:58:39 PM
> Subject: gromacs.org_gmx-users Digest, Vol 159, Issue 66
> 
> Send gromacs.org_gmx-users mailing list submissions to
>   gromacs.org_gmx-users@maillist.sys.kth.se
> 
> To subscribe or unsubscribe via the World Wide Web, visit
>   https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users
> or, via email, send a message with subject or body 'help' to
>   gromacs.org_gmx-users-requ...@maillist.sys.kth.se
> 
> You can reach the person managing the list at
>   gromacs.org_gmx-users-ow...@maillist.sys.kth.se
> 
> When replying, please edit your Subject line so it is more specific
> than "Re: Contents of gromacs.org_gmx-users digest..."
> 
> 
> Today's Topics:
> 
>1. Error in simulating graphene (B S Bhushan)
>2. Re: umbrella sampling (Justin Lemkul)
>3. Re: Concrete pull code explanation needed (Justin Lemkul)
>4. Re: ACETONITRILE with CHARM ff (Justin Lemkul)
>5. Re: Non-periodic COM Pulling (Justin Lemkul)
> 
> 
> --
> 
> Message: 1
> Date: Thu, 13 Jul 2017 15:33:16 +0530
> From: B S Bhushan 
> To: gromacs.org_gmx-users@maillist.sys.kth.se
> Subject: [gmx-users] Error in simulating graphene
> Message-ID:
>   
> Content-Type: text/plain; charset="UTF-8"
> 
> Dear Experts,
> 
> I am new to gromacs. I want to simulate electric double layer capacitance
> at graphene - liquid electrolyte interface. so far, I have practiced the
> tutorials from bevanlab page "Lysozyme in water" "biphasic systems". Next,
> I was trying to simulate the CNT/graphene tutorials available at thirdparty
> web pages. However I am getting errors while trying to generate topology
> file using *x2top* module.
> 
> When following the tutorial given at
> http://chembytes.wikidot.com/grocnt#toc
> I am getting the error,
> "*Inconsistency in user input:*
> *Could not find force field 'cnt_oplsaa' in current directory, install tree
> or*
> *GMXLIB path*."
> 
> 
> When following the tutorial given at
> http://machine-phase.blogspot.in/2009/04/single-wall-carbon-
> nanotubes-in-403.html
> I am getting the error,
> "*Fatal error:*
> *Could only find a forcefield type for 0 out of 168 atoms.*"
> 
> 
> Please suggest.
> Thank you so much for your time and knowledge.
> 
> 
> 
> 
> Awaiting suggestions,
> B S Bhushan
> Research Scholar,
> ABV - Indian Institute of Information Technology and Management, Gwalior,
> India.
> 
> 
> --
> 
> Message: 2
> Date: Thu, 13 Jul 2017 08:56:14 -0400
> From: Justin Lemkul 
> To: Discussion list for GROMACS users 
> Subject: Re: [gmx-users] umbrella sampling
> Message-ID: <995aa507-8875-8c41-674b-84eded63e...@vt.edu>
> Content-Type: text/plain; charset=utf-8; format=flowed

> > 
> > Dear Justin,
> > 
> > I increased the size of the simulation box to 4 nm. Indeed, the values of
> > pressure improved (averages around -1 bar or 1.5, 2 bar, or so). However,
> > the T keeps being overestimated (302 K). During NVT I got the right value,
> > so I decided to run the 200 ns NPT equilibration with Berendsen barostat
> > instead of P-R. I got the right T value, around 298.3 and a pressure of
> > 1.5 more or less. Then I launched a continuation test with 200 ps MD with
> > P-R (I couldn't use the -e option because it gives me the error Could not
> > find energy term named 'Box-Vel-XX', moreover, I didn't use P-R before so
> > I guess I shouldn't expect stored values for it?. But I am using -t ).
> > Despite I got low pressure averages around 0.5 bar, the temperature raised
> > again to 301-302. This happens very early in the simulation, which seems
> > to indicate that for sd integrator/thermostat and P-R barostat there is
> > something going on. I am getting practically the same density in between
> > simulation, so I guess I can say t
>  ha
> >   t in term of P, the system has been equilibrated. However, what can I do
> >   to get the right T for this system? If I was able to get the right T
> >   with Berendsen barostat, I don't understand what's wrong when I change
> >   to P-R.
> > 
> 
> Sounds buggy.  What version of GROMACS are you using?  There were temperature
> issues with the Langevin integrator in previous versions, but those should
> have
> been long since solved.

I am using version 5.1.2. I checked the release notes for version 5.1.3 and 
5.1.4 and the only thing related with sd integrator (none) or P-R barostat is 
this one 'Fixed Parrinello-Rahman with nstpcouple > 1'. Can this be the cause? 

Thank you.
> 
> -Justin
> 
-- 
Gromacs Users mailing list

* Please search the archive at 

Re: [gmx-users] Acetonitrile using CHARMM ff

[Justin Lemkul]

On 7/10/17 10:52 AM, Sonia Milena Aguilera Segura wrote:

Dear Justin,

Thank you for the answer. I changed the two parameters suggested in the mdp 
file and I ran again a minimization, 200 ps NVT, 200 ps NPT, and 1 ns MD for 
the two cases of the previous mail, and now I am getting densities around 771 
g/m3 which is slightly underestimated, but close to what other authors have 
obtained (774 others and 777 experimental). Also, I got slightly higher values 
for dielectric constants and diffusivities, also closer to another simulation 
with CHARMM. The energies also changed, but I guess that was expected. It looks 
like reproducing the dielectric constant with the current parameters is not 
possible. Is there anything that can be changed in order to get a better 
description? In terms of simulation, what is the dielectric constant depending 
of?



Structure, charge distribution, etc.  This may just be a limitation of an 
additive force field model.  We typically see neat liquid properties better 
reproduced with polarizable models.



Moreover, I observed that this time I got lower values for P during the NPT 
equilibration, but still is too far from 1 bar.  I really don't understand why 
for the NVT simulation I get a T around 298, but as soon as I turn on the 
pcoupl, the T rises to 300-301 K and the P gets average values of 7 and 4 bar 
(vs 8 and 14 for the previous simulations). Then at the end of the 1-ns MD the 
temperature remains around 301 and the P is -1 and 2.7 bar. Considering the 
parameters I am using, is there anything I can change to make the P coupling 
better? I am running a 3 nm box with 308 molecules. This is the full mdp file:



http://www.gromacs.org/Documentation/Terminology/Pressure

Your box is very small and will be subject to large fluctuations.

-Justin


; Run control
integrator   = sd   ; Langevin dynamics
tinit= 0
dt   = 0.0005
nsteps   = 200   ; 1 ns
nstcomm  = 100
;energygrps  = non-Water
; 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 = force-switch
rvdw-switch  = 1.0
rvdw = 1.2
; Apply long range dispersion corrections for Energy and Pressure
DispCorr  = no
; 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


Thank you very much,

Sonia Aguilera
PhD student
ENSCM

; Run control
integrator   = sd   ; Langevin dynamics
tinit= 0
dt   = 0.0005
nsteps   = 4000   ; 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


CHARMM uses a force switch, and only apply dispersion correction in the case of
lipid monolayers.

http://www.gromacs.org/Documentation/Terminology/Force_Fields/CHARMM

-Justin

==



--
==

Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 629
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalem...@outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul

Re: [gmx-users] Acetonitrile using CHARMM ff

[Sonia Milena Aguilera Segura]

Dear Justin, 

Thank you for the answer. I changed the two parameters suggested in the mdp 
file and I ran again a minimization, 200 ps NVT, 200 ps NPT, and 1 ns MD for 
the two cases of the previous mail, and now I am getting densities around 771 
g/m3 which is slightly underestimated, but close to what other authors have 
obtained (774 others and 777 experimental). Also, I got slightly higher values 
for dielectric constants and diffusivities, also closer to another simulation 
with CHARMM. The energies also changed, but I guess that was expected. It looks 
like reproducing the dielectric constant with the current parameters is not 
possible. Is there anything that can be changed in order to get a better 
description? In terms of simulation, what is the dielectric constant depending 
of?

Moreover, I observed that this time I got lower values for P during the NPT 
equilibration, but still is too far from 1 bar.  I really don't understand why 
for the NVT simulation I get a T around 298, but as soon as I turn on the 
pcoupl, the T rises to 300-301 K and the P gets average values of 7 and 4 bar 
(vs 8 and 14 for the previous simulations). Then at the end of the 1-ns MD the 
temperature remains around 301 and the P is -1 and 2.7 bar. Considering the 
parameters I am using, is there anything I can change to make the P coupling 
better? I am running a 3 nm box with 308 molecules. This is the full mdp file:

; Run control
integrator   = sd   ; Langevin dynamics
tinit= 0
dt   = 0.0005
nsteps   = 200   ; 1 ns
nstcomm  = 100
;energygrps  = non-Water
; 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 = force-switch
rvdw-switch  = 1.0
rvdw = 1.2
; Apply long range dispersion corrections for Energy and Pressure
DispCorr  = no
; 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


Thank you very much, 

Sonia Aguilera
PhD student
ENSCM
> ; Run control
> integrator   = sd   ; Langevin dynamics
> tinit= 0
> dt   = 0.0005
> nsteps   = 4000   ; 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

CHARMM uses a force switch, and only apply dispersion correction in the case of 
lipid monolayers.

http://www.gromacs.org/Documentation/Terminology/Force_Fields/CHARMM

-Justin

==
-- 
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Acetonitrile with CHARMM ff

[Justin Lemkul]

On 7/7/17 10:11 AM, Sonia Milena Aguilera Segura wrote:

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.302300.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. Constant20.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?  


I don't know the history of the acetonitrile parameters but very rarely to force 
field parameters reproduce every quantity of interest (heck, no one can even get 
water right).  There is a limit to how accurate a fixed-charge simulation can 
be, and perhaps you're finding that out.


Find the literature reference for the parameters and see what quantities are 
reproduced well and what the limitations might be.


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   = 4000   ; 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


CHARMM uses a force switch, and only apply dispersion correction in the case of 
lipid monolayers.


http://www.gromacs.org/Documentation/Terminology/Force_Fields/CHARMM

-Justin

--
==

Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 629
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalem...@outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul

==
--
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [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.302300.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. Constant20.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.565and 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   = 4000   ; 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=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 ]
> > ;  ijk  func   theta0   ktheta  ub0
> > kub
> > CG331CG1N1NG1T1 5   179.90   177.401600   0.
> > 0.00
> > 
> > [ dihedraltypes ]
> > ;  ijkl  func phi0 kphi  mult
> >   HGA3 CG331CG1N1  NG1T1 9 0.00 0.00 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
> 

-- 
Gromacs Users mailing list

* Please search the archive at 

Re: [gmx-users] Acetonitrile with CHARMM ff

[Justin Lemkul]

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=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 ]
;  ijk  func   theta0   ktheta  ub0 
 kub
CG331CG1N1NG1T1 5   179.90   177.401600   0. 
0.00

[ dihedraltypes ]
;  ijkl  func phi0 kphi  mult
  HGA3 CG331CG1N1  NG1T1 9 0.00 0.00 1
   
[ dihedraltypes ]

; 'improper' dihedrals
;  ijkl  func phi0 kphi
  HGA3 CG331CG1N1  NG1T1 9 0.00 0.00 1

Please let me know if this makes any sense. Is it correct to modify the 
angletype or should I only modify the initial structure?



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

--
==

Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 629
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalem...@outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul

==
--
Gromacs Users mailing list

* Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a 
mail to gmx-users-requ...@gromacs.org.

Re: [gmx-users] Acetonitrile with CHARMM ff

[Sonia Milena Aguilera Segura]

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=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 ]
;  ijk  func   theta0   ktheta  ub0 
 kub
CG331CG1N1NG1T1 5   179.90   177.401600   0. 
0.00

[ dihedraltypes ]
;  ijkl  func phi0 kphi  mult
 HGA3 CG331CG1N1  NG1T1 9 0.00 0.00 1 
  
[ dihedraltypes ]
; 'improper' dihedrals 
;  ijkl  func phi0 kphi
 HGA3 CG331CG1N1  NG1T1 9 0.00 0.00 1 

Please let me know if this makes any sense. Is it correct to modify the 
angletype or should I only modify the initial structure?

Is there someone who has already done simulations with acetonitrile? Any advice 
and feedback is very much welcome!!

Thank you very much, 

Have a nice day, 

Sonia Aguilera
PhD Student
ENSCM

> 
> Message: 2
> Date: Wed, 5 Jul 2017 17:47:08 -0400
> From: Justin Lemkul 
> To: gmx-us...@gromacs.org
> Subject: Re: [gmx-users] gromacs.org_gmx-users Digest, Acetonitrile
>   with CHARMM ff
> Message-ID: 
> Content-Type: text/plain; charset=utf-8; format=flowed
> 
> 
> 
> On 7/5/17 10:47 AM, Sonia Milena Aguilera Segura wrote:
> >> On 7/5/17 5:15 AM, Sonia Milena Aguilera Segura wrote:
> >>> Dear GROMACS users,
> >>>
> >>> I am currently interested in studying properties of some solvents, among
> >>> them acetonitrile and isopropanol. I would like to use CHARMM force field
> >>> for compatibility with other molecules and I am taking my initial
> >>> structure files from virtualchemistry.org. Does someone know how to run
> >>> the 6-sites model available with the CHARMM ff in gromacs? As I try to
> >>> run
> >>> the simulation box I get the error "No Defaul Proper Dih. Types". I
> >>> checked the ffbonded file and I didn't see the dihedrals defined as in
> >>> the
> >>> rtp file. And for what I've understood this is
> >>
> >> Dihedrals aren't required in .rtp files since pdb2gmx generates them.
> > 
> > 
> > Dear Justin,
> > 
> > Thank you for the comments. Yes, sorry, I was referring to the top file
> > generated by pdb2gmx instead of the rtp file. The top file has the
> > dihedrals defined whereas ffbonded files does not.
> > 
> > 
> >>
> >> because acetonitrile is a linear molecule and dihedrals for three colinear
> >> atoms
> >> this are mathematically
> >>
> >> There are still H-C-C-N dihedral terms.
> > 
> > 
> > Yes, they are named in the top file but not really defined in the ffbonded
> > file, which is why I get the error.
> > 
> 
> pdb2gmx will blindly generate all possible dihedrals; sometimes this isn't
> right.
> 
> It even appears that such parameters (H-C-C-N) do not exist in any CHARMM
> force
> field file.  Perhaps it's meant that way, but I have never tried to simulate
> anything with it.
> 
> > 
> >>
> >> undefined. Also, when I go to check the available itp for acetonitrile in
> >> virtualchemistry.org, I can see that there is a 7-sites model,with a dummy
> >> atom
> >> included. However, for this case, pdb and itp files do not match. I have
> >> seen
> >> that this can be sort of solved by fixing the angle as 179.9, but I really
> >> don't
> >> exactly what to change
> >>> or where in the force field files. I have no experience modifying
> >>> force
> >>> fields. Also, I've seen that for
> >>
> >> Don't modify the force field.  The CHARMM parameters for acetonitrile were
> >> generated in CHARMM, which can handle linear angles without the tricks
> >> that
> >> GROMACS requires with virtual sites.
> > 
> > 
> > So, what you are saying is that it is not possible to simulate the 6-sites
> > model of acetonitrile in GROMACS using CHARMM?
> 
> Never tried it.  AFAIK, GROMACS needs special treatment for linear species.
> Hopefully someone who has actually done such a simulation will say something,
> because my usefulness here is at its end :)
> 
> -Justin
> 
> > 
> >>
> >> the opls ff both pdb and itp files match, but I really need to use the
> >> CHARMM
> >> force field. Are the opls parameters compatible with CHARMM? Any advice?
> >>>
> >>
> >> Don't mix and match force fields.
> > 
> > 
> > Then, if I cannot run the 6-sites model, what's