[gmx-users] NEMD with sinusoidal pressure variation

[Aarushi Bhargava]

I am simulating a protein using gromacs 5.0.2 and want to apply time
varying pressure (preferably sine wave) on it. How should I do it in NEMD?
-- 
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] Drude force field

[Shraddha Parate]

Hi gmx users!

For polarizable effect, I was setting drude force field in Gromacs
I attempted running Lagrangian Dynamics simulation using Drude force field.
mdp file format was used from the supporting data in below paper.

[2015] Implementation of extended Lagrangian dynamics in GROMACS for
polarizable simulations using the classical Drude oscillator model.(
http://onlinelibrary.wiley.com/doi/10.1002/jcc.23937/abstract;jsessionid=E354008D7BE6AA1A02D26C061B49163A.f02t01
)




---

  ; RUN CONTROL PARAMETERS
  integrator = md-vv

  ; Start time and timestep in ps
  tinit = 0
  dt = 0.001
  nsteps = 1
  comm-mode = linear
  nstcomm = 100
  comm-grps = System

  ; OUTPUT CONTROL OPTIONS
  nstxout = 1000
  nstvout = 1000
  nstfout = 1000

  ; Output frequency for energies to log file and energy file
  nstlog = 100
  nstcalcenergy = 100
  nstenergy = 100

  ; NEIGHBORSEARCHING PARAMETERS
  cutoff-scheme = verlet
  nstlist = 10
  ns-type = Grid
  pbc = xyz
  rlist = 1.2

  ; OPTIONS FOR ELECTROSTATICS AND VDW
  coulombtype = PME
  rcoulomb = 1.2
  vdwtype = cutoff
  vdw-modifier = potential-switch
  rvdw-switch = 1.0
  rvdw = 1.2
  DispCorr = EnerPres

  ; TEMPERATURE COUPLING
  tcoupl = Nose-Hoover
  nsttcouple = 1
  nh-chain-length = 1

  ; Groups to couple separately
  tc-grps = Atoms Drues
  tau-t = 0.1 0.005
  ref-t = 298.15 1.0

  ; GENERATE VELOCITIES FOR STARTUP RUN, OTHERWISE DO
  ; NOT GENERATE NEW VELOCITIES IF PREVIOUSLY EQUILIBRATED
  gen-vel = yes
  gen-temp = 298.15
  gen-seed = -1

  ; OPTIONS FOR BONDS
  constraints = none ; can also be h-bonds, not all-bonds
  continuation = no


  ; DRUDE PARAMETERS
  drude = yes
  drude-mode = Lagrangian
  drude-t = 1.0
  drude-hardwall = yes
  drude-r = 0.02
  drude-tsteps = 5
---


I couldn't understand [tc-grps = Atoms Drude] option.
If the option is changed to [tc-grps = Protein SOL],I get a fatal error
message.

--
Fatal error:
Incorrect thermostat setup in relative_tstat, ti = -1
--

How to resolve this issue? Why the changed parameter gives a fatal error?

I look forward to your response.
Thank you for reading.
-- 
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] pdb2gmx

[‪farial tavakoli‬ ‪]

Hello
 I want to create a .gro and .top file from my protein that contains 379 
aminoacids in it's .pdb file by using gromos96 54a7 force field:
pdb2gmx -f protein.pdb -o protein.gro -water spce -ignh
but when gro and topology files are created , I see that message:
Start terminus GLY-12: GLY-NH3+
End terminus PRO-379: COO-
Checking for duplicate atoms
Generating any missing hydrogen atoms and/or adding termini.
Now there are 368 residues with 3856 atoms
Making bonds...
Number of bonds was 3938, now 3933
Generating angles, dihedrals and pairs...
Before cleaning: 6090 pairs
Before cleaning: 8285 dihedrals
Making cmap torsions...There are 2746 dihedrals, 2028 impropers, 5756 angles
  6090 pairs, 3933 bonds and 0 virtual sites
Total mass 42175.307 a.m.u.
Total charge -0.000 e
Writing topology

Back Off! I just backed up posre.itp to ./#posre.itp.46#

Writing coordinate file...

Back Off! I just backed up HDAC2.gro to ./#HDAC2.gro.1#
        - PLEASE NOTE 
You have successfully generated a topology from: HDAC2.pdb.
The Gromos54a7 force field and the spce water model are used.
        - ETON ESAELP 

However, i ignored this message and continued to use grompp
gmx grompp -f em.mdp -c solv.gro -p topol.top -o ions.tpr

but faced to this message:
  gmx grompp -f em.mdp -c solv.gro -p topol.top -o ions.tpr

Ignoring obsolete mdp entry 'title'

Back Off! I just backed up mdout.mdp to ./#mdout.mdp.29#
Setting the LD random seed to -1863558486
Generated 168 of the 1653 non-bonded parameter combinations
Excluding 3 bonded neighbours molecule type 'Protein_chain_A'
Excluding 2 bonded neighbours molecule type 'SOL'
Removing all charge groups because cutoff-scheme=Verlet
Analysing residue names:
There are:   368    Protein residues
There are: 11240  Water residues
Analysing Protein...
Number of degrees of freedom in T-Coupling group rest is 79005.00
Calculating fourier grid dimensions for X Y Z
Using a fourier grid of 72x72x72, spacing 0.114 0.114 0.114
Estimate for the relative computational load of the PME mesh part: 0.10
This run will generate roughly 3 Mb of data

Back Off! I just backed up ions.tpr to ./#ions.tpr.36#


Is there anyone to help me? Why pdb2gmx didnt consider all 379 amino acids of 
my protein and why GOMACS excluded 3 bonded neighours of the protein and 
SOL?and would you please tell what ' removing all charge groups because 
cutoff-scheme=Verlet ' means ? I t means the system is neutral ?

thanks in advanceFarial


-- 
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] Computing free energy differences from the dH/dl values directly

[Miguel Caro]

Hello,

I want to do thermodynamic integration directly from the ∂H/∂λ\partial H
/ \partial \lambda values in the md.xvg files, rather than using gmx
bar, because I want to do several manipulations (eventually). On the
tutorial on solvated methane how to use gmx bar is explained, however
this is a black box tool so I would like to understand a bit better what
it's doing and how it handles the raw derivatives.

I have a system (one methanol solvated in water) where I decouple
Coulomb interactions before decoupling vdW interactions. The md.xvg
files give me ∂H/∂λvdW\partial H / \partial \lambda_{vdW} and
∂H/∂λcou\partial H / \partial \lambda_{cou}, and so I am
computing ⟨∂H/∂λ⟩=⟨∂H/∂λvdW⟩∂λvdW/∂λ+⟨∂H/∂λc⟩∂λc/∂λ\langle \partial H /
\partial \lambda \rangle = \langle \partial H / \partial \lambda_{vdW}
\rangle \partial \lambda_{vdW} / \partial \lambda + \langle \partial H /
\partial \lambda_{c} \rangle \partial \lambda_{c} / \partial \lambda "by
hand".  I am then using a spline to generate smoother curves in λ\lambda
from the discrete array (21 data points, one for each λ\lambda value) of
the expectation values. I then integrate this smooth curve, although I
guess I could also use a quadrature rule applied to the original data
before smoothing. The final free energy difference I obtain like this is
about 5% larger than the value given by gmx bar.

I would like to know if what I am doing makes sense and if the
difference with gmx bar is to be expected (for instance because of how
the integral is performed or some other fundamental difference).

Many thanks,

Miguel

-- 

*Dr. Miguel Caro*
/Postdoctoral researcher/
Department of Electrical Engineering and Automation,
and COMP Centre of Excellence in Computational Nanoscience
Aalto University, Finland
Personal email: *mcar...@gmail.com*
Work: *miguel.c...@aalto.fi*
Website: http://mcaroba.dyndns.org
-- 
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] dihedral restraint

[Yu Ogawa]

Hi Gangotri,
You can increase the force constant.
Best, -Yu

2017-07-19 19:15 GMT+01:00 gangotri dey :

> Hello all,
>
> I am a relatively new user of gromacs. Hence, I am encountering some
> problem.
>
> I am studying few drug molecules onto a metal oxide surface. I would like
> to constrain the Ramachandran angles. I am using gromacs version 5.4.0.
>
> I have added the dihedral restraint angles in the itp file that reads as
> follows.
>
> [ dihedral_restraints ]
> ; ai   ajakal  type  phi  dphi  kfac
> 698 702 692 691 1 163.19 0 1
> 719 712 716 698 1 -84.24 0 1
>
> 712 716 698 689 1 -35.93 0 1
> 702 692 691 686 1 180.42 0 1
>
> I have read in this forum that there is no need to add any more key words
> like the previous version in the mdp file.
>
>
> However, the dihedral angles are changing and I am not sure how to restrain
> them.
> Is there any part that I am missing.
>
>
> *Thank you*
>
> *Gangotri Dey*
> --
> 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.
>
-- 
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] 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] pull COM "geometry = direction periodic" problem

[edesantis]

dear all,

I am studying the affinity between an antibody and an amyloid peptide;
I am interested in the evaluation of the PMF.

I tried to preform the COM pulling simulation using the "pull-geometry = 
direction", but for same the configurations (the ones in which the COMs 
of the two groups are quite far) I found the problem:
"Distance of pull group 1 is larger (...) than 0.49 times the box size ( 
... ). You might want to consider using "pull-geometry = 
direction-periodic" instead."


so according to the given suggestion, I tried to pull with the 
"pull-geometry = direction-periodic",

in this case I found the error :
"Fatal error:
184 particles communicated to PME rank 13 are more than 2/3 times the 
cut-off out of the domain decomposition cell of their charge group in 
dimension y.

This usually means that your system is not well equilibrated. "

the system is well equilibrated (according to the RMSD), the starting 
.gro file for gmx grompp is the same used for the geometry direction, 
and also the .mdp file is the same, except for the pull geometry 
directive,



so I don't understand what is the problem,


this is the COM part of the .mdp file used

; Pull code
pull= yes
pull_ngroups= 2
pull_ncoords= 1
pull_group1_name= Chain_Abeta
pull_group2_name= Chains_Antibody
pull_coord1_type= umbrella  ; harmonic biasing force
pull_coord1_geometry= direction-periodic
pull_coord1_groups  = 1 2
pull_coord1_vec =  38.207 68.611 29.8493
pull_coord1_rate= 0.002;
pull_coord1_k   = 1000  ; kJ mol^-1 nm^-2
pull_coord1_start   = yes   ; define initial COM distance > 
0


any of you have suggestions?


thanks in advance for your help
Emiliano



--
Emiliano De Santis
--
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.