Re: [gmx-users] Ethanol energies with CHARMM ff

[Sonia Milena Aguilera Segura]

>On 7/18/17 1:28 PM, Sonia Milena Aguilera Segura wrote:
>> Dear GROMACS users,
>> 
>> I am running a 4 nm box of ethanol and once I finished the MD, I got the 
>> right T, P, density. However, I noticed that my energies are odd. After 
>> several trials with different parameters and box sizes I >>am getting a 
>> Total Energy between -1.2 to 0.6 KJ/mol (already normalizaed, Total 
>> energy/#molecules). I checked the potential energy at the end of 
>> minimization and this was around -65 kJ/mol. Then, I >>observed that during 
>> NVT equilibration my potential energy decreased to a value around -33 kJ/mol 
>> with a kinetic energy also around the same value (33 kJ/mol), and that's why 
>> my final total energies >>and, therefore, enthalpies are giving values 
>> around 0. I ran 

>I don't understand how you're calculating your enthalpies.  The kinetic energy 
>is irrelevant because it should cancel between the gas and liquid phases for a 
>given temperature, so the dHvap value is /N +  + RT

>>also simulations with water, isopropanol, and acetonitrile and I am getting 
>>values of -52, -206, and-48, which seem reasonable to me. My paremeters have 
>>been already discussed in a previous mail (please see acetonitrile

>Reasonable based on...?
>
>-Justin

Dear Justin,

Currently I am not trying to calculate dHvap, only comparing the energies 
(Total energy/#molecules and enthalphy with -nmol and -fluct_props option for 
g_energy, that gives me enthalpies in kJ/mol). If I am not mistaken, in 
practice this enthalpy is calculated as the total energy (internal energy) plus 
de pV correction, all divided by the number of molecules in the system to give 
a value of enthalpy in kJ/mol. Is there any way to compare this energies 
against experimental values? For what I've understood, this value can be only 
comparable against a computational chemistry calculation (DFT, etc) to compare 
the total energy of the system. Please correct me if I am wrong. 

When I say that the values I am getting seem reasonable I meant that they are 
all non-zero values, and moreover, without really know if the value I am 
getting correspond to any experimental value and/or ab initio calculation, at 
least they are all negative and they seem to show and correspond to the 
strength of the molecular interactions, e.g. isopropanol displays higher 
energies than water and acetonitrile, which kind of correspond to strength of 
their molecular interactions. Moreover, the problem comes when I see the values 
of ethanol. No matter what, this value shouldn't be zero. Or am I missing 
something? 

Perhaps I should have started by asking the physical meaning of this values and 
the right approach to analyze the energies of the system. Is comparing the 
energies as described before a bad practice? Is computing the dHvap the only 
and right approach?.  In this case I would also need to simulate the gas phase 
for my solvents, right?. However, all of them are liquid at 298. How can I 
simulate the gas phase at this T?

Thank you very much in advance, 

Sonia Aguilera
PhD student
ENSCM



>>with CHARMM ff), and they seem to be right to be used with CHARMM ff. I am 
>>using 
>>version 2016.3 (I did the sa
>>   me in 5.1.2), sd integrator, Berendsen barostat for equilibration and P-R 
>> for MD production. Any ideas of what should I look for or what can I be 
>> doing wrong?
> 
-- 
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.

[gmx-users] Ethanol energies with CHARMM ff

[Sonia Milena Aguilera Segura]

Dear GROMACS users, 

I am running a 4 nm box of ethanol and once I finished the MD, I got the right 
T, P, density. However, I noticed that my energies are odd. After several 
trials with different parameters and box sizes I am getting a Total Energy 
between -1.2 to 0.6 KJ/mol (already normalizaed, Total energy/#molecules). I 
checked the potential energy at the end of minimization and this was around -65 
kJ/mol. Then, I observed that during NVT equilibration my potential energy 
decreased to a value around -33 kJ/mol with a kinetic energy also around the 
same value (33 kJ/mol), and that's why my final total energies and, therefore, 
enthalpies are giving values around 0. I ran also simulations with water, 
isopropanol, and acetonitrile and I am getting values of -52, -206, and-48, 
which seem reasonable to me. My paremeters have been already discussed in a 
previous mail (please see acetonitrile with CHARMM ff), and they seem to be 
right to be used with CHARMM ff. I am using version 2016.3 (I did the sa
 me in 5.1.2), sd integrator, Berendsen barostat for equilibration and P-R for 
MD production. Any ideas of what should I look for or what can I be doing 
wrong? 

Thank you very much!!!, 

Sonia Aguilera 
PhD student 
ENSCM 
-- 
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] sd integrator and P-R barostat/ACETONITRILE with CHARMM ff

[Sonia Milena Aguilera Segura]

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

Dear Justin, 

I upgraded to 2016.3 and did several tests. In all cases I ran first a NVT with 
sd and 200 ps NPT with Berendsen. I always got a T of 297.6-298.6 for 
equilibrations. Then, for instance,  I ran a 4 nm box with isopropanol, which 
is bigger than acetonitrile, and got a T of 300 K instead of 302 K, for the 
same system with the 5.1.2 version. For a 4 nm box of acetonitrile I got also a 
small improvement: 299.5 vs 300-301 for the previous version. So, I guess there 
was something wrong with the version. Still, I am wondering if the level of 
precision is enough (an error of +-1.5 K), considering that I am getting very 
precise values for Berendsen barostat, even for the 3 nm boxes for both 
solvents. Also, my average P is -1.5 and 1.5.. which is giving me very similar 
results for density in all cases. Can I assume that those values are ok (even 
the negative) considering the expected fluctuations during the MD simulation?

Thank you, 

Best regards, 

Sonia Aguilera 
PhD student
ENSCM
-- 
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 with CHARM ff

[Sonia Milena Aguilera Segura]

> 
> Message: 3
> Date: Tue, 11 Jul 2017 12:57:05 -0400
> From: Justin Lemkul 
> > 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.

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

Thank you in advance, 

Sonia Aguilera 
PhD student
ENSCM

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

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

[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

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 <jalem...@vt.edu>
> To: gmx-us...@gromacs.org
> Subject: Re: [gmx-users] gromacs.org_gmx-users Digest, Acetonitrile
>   with CHARMM ff
> Message-ID: <c935a574-ffc8-cf5f-3a3c-3131d906a...@vt.edu>
> 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

Re: [gmx-users] gromacs.org_gmx-users Digest, Acetonitrile with CHARMM ff

[Sonia Milena Aguilera Segura]

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


> 
> 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? 

> 
> 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 the advice? I would really need 
to run acetonitrile with CHARMM, because I have other series of molecules such 
as cellulose and other saccharides. In that case, what would be the most 
compatible for field to simulate organic solvents and saccharides?  

> 
> > Also, I cannot find an already optimized structure with all hydrogens for
> > isopropanol. Could somenone provide one?
> 
> Take any valine side chain and use it.  You can generate any missing H easily
> with a suitable .hdb entry.
> 
> -Justin

Thank you for your advice, 

Sonia Aguilera
PhD student
ENSCM
-- 
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.

[gmx-users] Acetonitrile and isopropanol with CHARMM ff

[Sonia Milena Aguilera Segura]

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 because acetonitrile is a linear molecule and dihedrals for 
three colinear atoms this are mathematically 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 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? 

Also, I cannot find an already optimized structure with all hydrogens for 
isopropanol. Could somenone provide one? 

Thank you very much in advance for your time and help!! 

Best regards, 

Sonia Aguilera 
PhD student 
ENSCM 
Montpellier, France 

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