Re: [gmx-users] regarding tha langevin thermostat (integrator sd)

[Mark Abraham]

Hi,

On Thu, May 11, 2017 at 11:45 PM Felix Y Yang  wrote:

> Hello,
>
> I am trying to use a Langevin Thermostat to only heat up a specific group
> of atoms in my molecule specified by an index file, while not doing
> anything to the rest of the system (grouped as another group in the index
> file). Is there any way I can achieve this effect?
>

Either the sd integrator or the v-rescale stochastic thermostat support
temperature coupling groups, which in effect functions as a langevin
thermostat. You can disable coupling to any group by choosing the
appropriate time constant...
http://manual.gromacs.org/documentation/2016.3/user-guide/mdp-options.html#temperature-coupling

Can I set up temperature annealing and set ref_t for the rest of the
> system as constant 300K over the run period (which is the temperature
> achieved for the system prior via v-rescale thermostat), while heating up
> the specific group to a higher temperature over the same time? Would this
> have the same effect as not affecting the rest of my system?
>

Yes, that kind of approach will be useful if you want gradual heating. But
I'm not sure whether you're correctly describing "not affecting the rest of
the system" if you actually mean leaving a normal coupling group on it ;-)

Mark


> Any feedback would be greatly appreciated. Thanks.
>
> Felix Yang
> Graduate Master's Student
> Department of Chemistry & Biochemistry
> University of California, San Diego
>
> --
> 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.

[gmx-users] regarding tha langevin thermostat (integrator sd)

[Felix Y Yang]

Hello,

I am trying to use a Langevin Thermostat to only heat up a specific group
of atoms in my molecule specified by an index file, while not doing
anything to the rest of the system (grouped as another group in the index
file). Is there any way I can achieve this effect?

Can I set up temperature annealing and set ref_t for the rest of the
system as constant 300K over the run period (which is the temperature
achieved for the system prior via v-rescale thermostat), while heating up
the specific group to a higher temperature over the same time? Would this
have the same effect as not affecting the rest of my system?

Any feedback would be greatly appreciated. Thanks.

Felix Yang
Graduate Master's Student
Department of Chemistry & Biochemistry
University of California, San Diego

-- 
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] RMSF problem

[Franco Henrique]

Dear Justin Lemkul,

Thank’s for prompt reply.

My protein has 198 aminoacides (PDB: 2GOJ; 
http://www.rcsb.org/pdb/explore.do?structureId=2GOJ). After RMSF command (gmx 
rmsf -f .xtc -s .tpr -o .xvg -res), the output has only 103 aminoacides as 
visualized by Grace software. I am running gromacs 5.1.2 in ubuntu.

Any inquiries are more than welcome.

Best regards,
Franco Henrique.
--
Doutor em Biotecnologia
Professor Assistente B 
Departamento de Saúde - DSAU
Universidade Estadual de Feira de Santana - UEFS
Avenida Transnordestina s/n, bairro Novo Horizonte CEP 44.036-900 - Feira de 
Santana – BA

Em 11 de mai de 2017, à(s) 14:28, Justin Lemkul  escreveu:



On 5/11/17 1:21 PM, Franco Henrique wrote:
> Dear,
> 
> I am have problem with RMSF data, because after I apply this command-line: 
> gmx rmsf -f xxx.xtc -s yyy.tpr -o zzz.xvg -res, the .xvg file has not all 
> residues of protein. If possible, I would like to know a tip to solve this 
> problem.
> 

How many residues do you have, and how many did the program output?  No GROMACS 
program will magically get rid of anything; it will analyze what you tell it 
to, so without more information about what you did, which group you chose, etc. 
it's hard to suggest anything.

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

-- 
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] RMSF problem

[Justin Lemkul]

On 5/11/17 1:21 PM, Franco Henrique wrote:

Dear,

I am have problem with RMSF data, because after I apply this command-line: gmx 
rmsf -f xxx.xtc -s yyy.tpr -o zzz.xvg -res, the .xvg file has not all residues 
of protein. If possible, I would like to know a tip to solve this problem.



How many residues do you have, and how many did the program output?  No GROMACS 
program will magically get rid of anything; it will analyze what you tell it to, 
so without more information about what you did, which group you chose, etc. it's 
hard to suggest anything.


-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] ligand moving out during umbrella sampling

[Justin Lemkul]

On 5/11/17 1:25 PM, abhisek Mondal wrote:

Alright. I'm trying as per your advice. Actually my system is pretty big
and due to constraint of computation power it is taking me more time to
test vector setup. I really appreciate your time. Thanks a lot for your
suggestions.

One thing I'm worried about here. In previous mail, I mentioned the
"distance at start" and "ref at t=0" is negative. What does the negative
value signify, I mean how can distance be negative ? Is is due to my poorly
given vector definition or something else ?



It's a vector quantity.  It can be positive or negative depending on how you 
define the groups and what their positions are.  And you're telling grompp to 
pull along the negative-z dimension.


-Justin


On Thu, May 11, 2017 at 9:03 PM, Justin Lemkul  wrote:




On 5/11/17 9:21 AM, abhisek Mondal wrote:


On Thu, May 11, 2017 at 11:02 AM, abhisek Mondal 
wrote:

Hi,


Thank you for the explanation. It really cleared some concepts. But I'm
still having my ligand moving in this step. I have modified the code as:
; Pull code
pull= umbrella
pull_ngroups= 1
pull_group0 = Protein_chain_A
pull_group1 = ACO
pull_geometry   = direction ; simple distance increase
pull_dim   = Y Y Y ; not to allow ligand move along other
dir
pull_rate1 = 0.0
pull_k1   = 1000   ; kJ mol^-1 nm^-2
pull_start   = yes   ; define initial COM distance > 0
pull_vec1   = 0 0 -1



Note that with "direction" geometry, only pull_vec1 is acting.  pull_dim
is ignored.

The ligand was previously moving along x,y direction when I was using

pull_dim  = N N Y. So I changed it to Y in all direction and provided 0
as
vector  and pull_rate1=0.0, so that it does not move much. But at the end
of a 10ns run, I see that the ligand is still moving as it was earlier.

It shows me:

 Pull group  natoms  pbc atom  distance at start reference at t=0
   0  1132936665
   159  1618  -1.555-1.555
Is it ok withe negative value ? Anyway this setup is not working.



Again you're trying to just apply the restraint to one dimension and it
looks to be fairly arbitrary.  I already suggested using the vector
connecting the ligand COM with the binding site residues' COM and using
that as pull_vec1.  Draw it out.  It makes a lot more sense than trying to
restrain only along one axis, which as I have said before, makes no sense
in this case.

-Justin



I have choose my reaction coordinate to be along -Z axis and want to
apply
biasing potential accordingly with restraining the ligand movement. Can
you
please suggest where am I failing with this code ?

Thank you.



On Tue, May 9, 2017 at 1:11 AM, Justin Lemkul  wrote:




On 5/8/17 10:00 AM, abhisek Mondal wrote:

On Sun, May 7, 2017 at 11:37 PM, Justin Lemkul  wrote:





On 5/7/17 1:57 AM, abhisek Mondal wrote:

Hi,



For your ease of understanding regarding what is happening during
this
above said umbrella-mdrun, I have shared the trajectory video file
the
following link.
https://drive.google.com/drive/folders/0B6O-L5Y7BiGJQ1FIc2tIRFE2dE0

Is this normal given that the mdp code being used ? I basically have
no
idea with this step, so please help me out. I'm using gromacs-4.6.2.


Your setup is incorrect.  You're applying a biasing potential only


along
z, so the ligand can move freely along x and y.  A protein-ligand
complex
has spherical symmetry, so you should set the reaction coordinate to
the
vector connecting the ligand with some suitable subset of interacting
protein residues.




I don't get it.
We are trying not to move our configuration (generated after pulling
simulation) along the reaction coordinate, so for restraining we are
supposed to set pull_rate1=0.0.



Of course.  But you said set pull_k = 0 which does not make sense.  The
pulling rate *is* zero during umbrella sampling (no net displacement,
restrain to the specified distance along the reaction coordinate) and
pulling force constant should be non-zero.

If applying biasing potential only along z is causing movement along x
and


y then what if we apply the biasing potential along x,y,z ? Will it
cause
any good in restraining the ligand?


This is how it should be done.  The reaction coordinate should be

suitably defined based on the geometry of the system.  As I suggested
before, choose some representative residues in the active site as one
group
and the ligand as the other.  Thus defines the reaction coordinate
without
any presupposition of anything being aligned with a Cartesian axis,
which
is rarely the case.

Moreover, you said previously "A protein-ligand complex has spherical


symmetry, so you should set the reaction coordinate to the vector
connecting the ligand with some suitable subset of interacting protein
residues.". It is really unclear to me, could you 

Re: [gmx-users] ligand moving out during umbrella sampling

[abhisek Mondal]

Alright. I'm trying as per your advice. Actually my system is pretty big
and due to constraint of computation power it is taking me more time to
test vector setup. I really appreciate your time. Thanks a lot for your
suggestions.

One thing I'm worried about here. In previous mail, I mentioned the
"distance at start" and "ref at t=0" is negative. What does the negative
value signify, I mean how can distance be negative ? Is is due to my poorly
given vector definition or something else ?

On Thu, May 11, 2017 at 9:03 PM, Justin Lemkul  wrote:

>
>
> On 5/11/17 9:21 AM, abhisek Mondal wrote:
>
>> On Thu, May 11, 2017 at 11:02 AM, abhisek Mondal 
>> wrote:
>>
>> Hi,
>>>
>>> Thank you for the explanation. It really cleared some concepts. But I'm
>>> still having my ligand moving in this step. I have modified the code as:
>>> ; Pull code
>>> pull= umbrella
>>> pull_ngroups= 1
>>> pull_group0 = Protein_chain_A
>>> pull_group1 = ACO
>>> pull_geometry   = direction ; simple distance increase
>>> pull_dim   = Y Y Y ; not to allow ligand move along other
>>> dir
>>> pull_rate1 = 0.0
>>> pull_k1   = 1000   ; kJ mol^-1 nm^-2
>>> pull_start   = yes   ; define initial COM distance > 0
>>> pull_vec1   = 0 0 -1
>>>
>>>
> Note that with "direction" geometry, only pull_vec1 is acting.  pull_dim
> is ignored.
>
> The ligand was previously moving along x,y direction when I was using
>>> pull_dim  = N N Y. So I changed it to Y in all direction and provided 0
>>> as
>>> vector  and pull_rate1=0.0, so that it does not move much. But at the end
>>> of a 10ns run, I see that the ligand is still moving as it was earlier.
>>>
>>> It shows me:
>>  Pull group  natoms  pbc atom  distance at start reference at t=0
>>0  1132936665
>>159  1618  -1.555-1.555
>> Is it ok withe negative value ? Anyway this setup is not working.
>>
>>
> Again you're trying to just apply the restraint to one dimension and it
> looks to be fairly arbitrary.  I already suggested using the vector
> connecting the ligand COM with the binding site residues' COM and using
> that as pull_vec1.  Draw it out.  It makes a lot more sense than trying to
> restrain only along one axis, which as I have said before, makes no sense
> in this case.
>
> -Justin
>
>
>>> I have choose my reaction coordinate to be along -Z axis and want to
>>> apply
>>> biasing potential accordingly with restraining the ligand movement. Can
>>> you
>>> please suggest where am I failing with this code ?
>>>
>>> Thank you.
>>>
>>>
>>>
>>> On Tue, May 9, 2017 at 1:11 AM, Justin Lemkul  wrote:
>>>
>>>

 On 5/8/17 10:00 AM, abhisek Mondal wrote:

 On Sun, May 7, 2017 at 11:37 PM, Justin Lemkul  wrote:
>
>
>
>> On 5/7/17 1:57 AM, abhisek Mondal wrote:
>>
>> Hi,
>>
>>>
>>> For your ease of understanding regarding what is happening during
>>> this
>>> above said umbrella-mdrun, I have shared the trajectory video file
>>> the
>>> following link.
>>> https://drive.google.com/drive/folders/0B6O-L5Y7BiGJQ1FIc2tIRFE2dE0
>>>
>>> Is this normal given that the mdp code being used ? I basically have
>>> no
>>> idea with this step, so please help me out. I'm using gromacs-4.6.2.
>>>
>>>
>>> Your setup is incorrect.  You're applying a biasing potential only
>>>
>> along
>> z, so the ligand can move freely along x and y.  A protein-ligand
>> complex
>> has spherical symmetry, so you should set the reaction coordinate to
>> the
>> vector connecting the ligand with some suitable subset of interacting
>> protein residues.
>>
>>
>
> I don't get it.
> We are trying not to move our configuration (generated after pulling
> simulation) along the reaction coordinate, so for restraining we are
> supposed to set pull_rate1=0.0.
>
>
 Of course.  But you said set pull_k = 0 which does not make sense.  The
 pulling rate *is* zero during umbrella sampling (no net displacement,
 restrain to the specified distance along the reaction coordinate) and
 pulling force constant should be non-zero.

 If applying biasing potential only along z is causing movement along x
 and

> y then what if we apply the biasing potential along x,y,z ? Will it
> cause
> any good in restraining the ligand?
>
>
> This is how it should be done.  The reaction coordinate should be
 suitably defined based on the geometry of the system.  As I suggested
 before, choose some representative residues in the active site as one
 group
 and the ligand as the other.  Thus defines the reaction coordinate
 without
 any presupposition of anything being aligned with a 

[gmx-users] RMSF problem

[Franco Henrique]

Dear,

I am have problem with RMSF data, because after I apply this command-line: gmx 
rmsf -f xxx.xtc -s yyy.tpr -o zzz.xvg -res, the .xvg file has not all residues 
of protein. If possible, I would like to know a tip to solve this problem.

Best regards,
Franco Henrique.
--
Doutor em Biotecnologia
Professor Assistente B 
Departamento de Saúde - DSAU
Universidade Estadual de Feira de Santana - UEFS
Avenida Transnordestina s/n, bairro Novo Horizonte CEP 44.036-900 - Feira de 
Santana – BA

-- 
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] Regarding simulation radioactive material

[Justin Lemkul]

Again, please keep the discussion on the mailing list.

On 5/10/17 9:47 PM, RAHUL SURESH wrote:

Dear Justin

Thank you. I have to think over it a lot. To reparametarise a whole protein, I
do not think it's possible or an easy job. The next pull down is the Change is
mass should correspond to change in neutrons and should affect the stability of
an atom which I don't whether it will work in any right way.

Indeed, please let me know any possible way to reperameterize a molecule. I can
simply try it with some single amino acids.



The way to parametrize anything is to follow the same protocol as the force 
field you're trying to modify.  In your case, the bonded terms would be 
essential because vibrational frequencies will change.  I doubt there will be 
differences in the nonbonded terms and doing anything like this strikes me as 
wasted effort because the MM approximation will likely view C-12 and C-14 as 
basically the same from the perspective of nonbonded interactions.


-Justin



On Thu, 11 May 2017 at 1:35 A
M, Justin Lemkul > wrote:


Please keep the discussion on the gmx-users mailing list.  I am not a 
private
help service.

On 5/10/17 2:58 PM, RAHUL SURESH wrote:
> Dear Justin
>
> I would like to simulate a radioactive peptide.. I think this is the  
first of
> its kind.  I wish to study interaction between Radioactive peptide and a
protein.
>
> To make it simple if I want to convert all my carbon into radioactive 
isotopes
> c14, and simulate is it possible ?
>

Possible, yes, but not correct.  If you change the mass of any atoms, you 
change
their vibrational frequencies, which means you need to reparametrize all the
bonded terms involving the modified atom types.  That means bonds, angles, 
and
dihedrals.  Whether or not this gives anything meaningful in terms of 
modeling
the interactions of a radioactive peptide with a larger protein is another
question.  If you're not changing any elements of the nonbonded terms 
(which you
likely shouldn't) then I doubt the interactions will be different in any
statistically meaningful way from the default parameters.  Re-doing bonded
parameters is a LOT of work.

-Justin

> References, I have found few  experimental papers on this. Where they have
Used
> isotopes of hydrogen. Tritium instead of hydrogen.
>
> If any information regarding this is available, will be a great use for 
me.
>
> Thank you.
> Sorry for writing in person
> --
> */Regards,/*
> */Rahul Suresh/*
> */Research Scholar/*
> */Bharathiar University/*
> */Coimbatore/*

--
==

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

==

--
*/Regards,/*
*/Rahul Suresh/*
*/Research Scholar/*
*/Bharathiar University/*
*/Coimbatore/*


--
==

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] Using the md integrator for calculating free energy of solvation

[Justin Lemkul]

On 5/11/17 9:49 AM, Dan Gil wrote:

Thank you Dr. Lemkul,

The simulations ran, but I have some a question about the results.

Is it natural to have a nonzero mass_lambda dH/dL? My intuition was that
the potential energy does not depend on mass, and my thinking was that the
mass contribution should be very small, if not zero.



It should.  Maybe it's a convergence issue?

-Justin


Best Regards,

Dan

On Wed, May 3, 2017 at 10:46 AM, Justin Lemkul  wrote:




On 5/3/17 10:43 AM, Dan Gil wrote:


Hi Dr. Kausar,

If I want the molecule to have all of its nonbonded interactions with the
rest of the system (solvent) at each lambda, don't I want to have vdwq
option ON for both couple-lambda0 and couple-lambda1?

What I am attempting to do is mutate molecule A to molecule B through the
simulation.



Relative free energy calculations are a little bit of a strange beast; you
do need to set the B-state (couple-lambda1) to "none" and the B-state types
and charges are read from [atoms].


-Justin

Best Regards,


Dan

On Wed, May 3, 2017 at 1:04 AM, Tasneem Kausar  wrote:

Thank you Dr. Lemkul,


Following your advice, I was able to calculate the solvation free energy


of


two molecules that I am interested in. The results, unfortunately, were


not


what I expected. I want to try an alternative method, that is, the
alchemical thermodynamic integration method.

If I understand this correctly, I should be able to calculate the free
energy difference between two systems A and B as I change the topology
of
the molecule from A to B in small steps.

I wasn't able to find a tutorial online, but I attempted to try the


method


and I obtained this error from grompp:

WARNING 1 [file grompp.mdp, line 43]:
  The lambda=0 and lambda=1 states for coupling are identical

I was hoping if you could help me understand what I am missing in my
topology or mdp file.

Topology:
 [ moleculetype ]
HEPT 3

 [ atoms ]
 1   opls_135  1   HEPTCH3  1 -0.180 12.011
 opls_9610.36012.011
 2   opls_136  1   HEPTCH2  2 -0.120 12.011
 opls_9620.24012.011
 3   opls_136  1   HEPTCH2  3 -0.120 12.011
 opls_9620.24012.011
 4   opls_136  1   HEPTCH2  4 -0.120 12.011
 opls_9620.24012.011
 5   opls_136  1   HEPTCH2  5 -0.120 12.011
 opls_9620.24012.011
 6   opls_136  1   HEPTCH2  6 -0.120 12.011
 opls_9620.24012.011
 7   opls_135  1   HEPTCH3  7 -0.180 12.011
 opls_9610.36012.011
 8   opls_140  1   HEPT  H  1  0.060  1.008
 opls_965   -0.12018.998
 9   opls_140  1   HEPT  H  1  0.060  1.008
 opls_965   -0.12018.998
10   opls_140  1   HEPT  H  1  0.060  1.008
 opls_965   -0.12018.998
11   opls_140  1   HEPT  H  2  0.060  1.008
 opls_965   -0.12018.998
12   opls_140  1   HEPT  H  2  0.060  1.008
 opls_965   -0.12018.998
13   opls_140  1   HEPT  H  3  0.060  1.008
 opls_965   -0.12018.998
14   opls_140  1   HEPT  H  3  0.060  1.008
 opls_965   -0.12018.998
15   opls_140  1   HEPT  H  4  0.060  1.008
 opls_965   -0.12018.998
16   opls_140  1   HEPT  H  4  0.060  1.008
 opls_965   -0.12018.998
17   opls_140  1   HEPT  H  5  0.060  1.008
 opls_965   -0.12018.998
18   opls_140  1   HEPT  H  5  0.060  1.008
 opls_965   -0.12018.998
19   opls_140  1   HEPT  H  6  0.060  1.008
 opls_965   -0.12018.998
20   opls_140  1   HEPT  H  6  0.060  1.008
 opls_965   -0.12018.998
21   opls_140  1   HEPT  H  7  0.060  1.008
 opls_965   -0.12018.998
22   opls_140  1   HEPT  H  7  0.060  1.008
 opls_965   -0.12018.998
23   opls_140  1   HEPT  H  7  0.060  1.008
 opls_965   -0.12018.998
... And so on with the bonds, pairs, angles, and dihedrals directive.

MDP File:

;Integration Method and Parameters
integrator   = md
nsteps   = 1
dt = 0.002
nstenergy= 1000
nstlog   = 5000

;Output Control
nstxout = 0
nstvout = 0

;Cutoff Schemes
cutoff-scheme= 

Re: [gmx-users] ligand moving out during umbrella sampling

[Justin Lemkul]

On 5/11/17 9:21 AM, abhisek Mondal wrote:

On Thu, May 11, 2017 at 11:02 AM, abhisek Mondal 
wrote:


Hi,

Thank you for the explanation. It really cleared some concepts. But I'm
still having my ligand moving in this step. I have modified the code as:
; Pull code
pull= umbrella
pull_ngroups= 1
pull_group0 = Protein_chain_A
pull_group1 = ACO
pull_geometry   = direction ; simple distance increase
pull_dim   = Y Y Y ; not to allow ligand move along other
dir
pull_rate1 = 0.0
pull_k1   = 1000   ; kJ mol^-1 nm^-2
pull_start   = yes   ; define initial COM distance > 0
pull_vec1   = 0 0 -1



Note that with "direction" geometry, only pull_vec1 is acting.  pull_dim is 
ignored.


The ligand was previously moving along x,y direction when I was using
pull_dim  = N N Y. So I changed it to Y in all direction and provided 0 as
vector  and pull_rate1=0.0, so that it does not move much. But at the end
of a 10ns run, I see that the ligand is still moving as it was earlier.


It shows me:
 Pull group  natoms  pbc atom  distance at start reference at t=0
   0  1132936665
   159  1618  -1.555-1.555
Is it ok withe negative value ? Anyway this setup is not working.



Again you're trying to just apply the restraint to one dimension and it looks to 
be fairly arbitrary.  I already suggested using the vector connecting the ligand 
COM with the binding site residues' COM and using that as pull_vec1.  Draw it 
out.  It makes a lot more sense than trying to restrain only along one axis, 
which as I have said before, makes no sense in this case.


-Justin



I have choose my reaction coordinate to be along -Z axis and want to apply
biasing potential accordingly with restraining the ligand movement. Can you
please suggest where am I failing with this code ?

Thank you.



On Tue, May 9, 2017 at 1:11 AM, Justin Lemkul  wrote:




On 5/8/17 10:00 AM, abhisek Mondal wrote:


On Sun, May 7, 2017 at 11:37 PM, Justin Lemkul  wrote:




On 5/7/17 1:57 AM, abhisek Mondal wrote:

Hi,


For your ease of understanding regarding what is happening during this
above said umbrella-mdrun, I have shared the trajectory video file the
following link.
https://drive.google.com/drive/folders/0B6O-L5Y7BiGJQ1FIc2tIRFE2dE0

Is this normal given that the mdp code being used ? I basically have no
idea with this step, so please help me out. I'm using gromacs-4.6.2.


Your setup is incorrect.  You're applying a biasing potential only

along
z, so the ligand can move freely along x and y.  A protein-ligand
complex
has spherical symmetry, so you should set the reaction coordinate to the
vector connecting the ligand with some suitable subset of interacting
protein residues.




I don't get it.
We are trying not to move our configuration (generated after pulling
simulation) along the reaction coordinate, so for restraining we are
supposed to set pull_rate1=0.0.



Of course.  But you said set pull_k = 0 which does not make sense.  The
pulling rate *is* zero during umbrella sampling (no net displacement,
restrain to the specified distance along the reaction coordinate) and
pulling force constant should be non-zero.

If applying biasing potential only along z is causing movement along x and

y then what if we apply the biasing potential along x,y,z ? Will it cause
any good in restraining the ligand?



This is how it should be done.  The reaction coordinate should be
suitably defined based on the geometry of the system.  As I suggested
before, choose some representative residues in the active site as one group
and the ligand as the other.  Thus defines the reaction coordinate without
any presupposition of anything being aligned with a Cartesian axis, which
is rarely the case.

Moreover, you said previously "A protein-ligand complex has spherical

symmetry, so you should set the reaction coordinate to the vector
connecting the ligand with some suitable subset of interacting protein
residues.". It is really unclear to me, could you please give me some
examples to understand it more simply? I had pulled the ligand along -z
axis, doesn't it mean that the reaction coordinate is to be that way ?
The
fact that I'm struggling with is to restrain the pull configurations for
further sampling.



The reaction coordinate is whatever you define it to be.  Whether or not
pulling along the z-axis makes sense depends on the orientation of the
system and the intrinsic geometry.  In your case, it doesn't make sense.
In my case (the tutorial, the unidirectional growth of an amyloid fibril)
it does make sense to use a single Cartesian axis for the SMD portion and
subsequent umbrella sampling.

I'm really a beginner, so maybe I'm asking stupid questions. Please give

me
some advise. I'm really unable to decipher the scenario in comparison to
your amyloid 

Re: [gmx-users] Inconsistence of SPC/E and TIP4P water in Gromacs

[Justin Lemkul]

On 5/11/17 1:53 AM, Mark Abraham wrote:

Hi,

On Thu, May 11, 2017 at 3:31 AM ZUO Taisen  wrote:



Hi guys:

 I have compared the SPC/E and TIP4P water of the opls force field in
Gromacs.But there is something inconsistent of the potential energy in the
system


The potential energy of SPC/E water is -4.71400e+04/1000=-47.14kJ/mol
which is far from the results from literature (Peter G. Kusalik,
Science,vol,265,26,p1219-1221,1994)
The potential energy TIP4P water is -4.17875e+04/1000=41.787kJ/mol which
is very close to the literature -41.8kJ/mol(Peter G. Kusalik,
Science,vol,265,26,p1219-1221,1994)


I have attached all the files simulating the SPC/E and TIP4P water(.trr
.xtc files are delected) and also the related paper in my last email but
was suspended because of the files were too big.So no files was attached
this time but the .mdp file was posted as following.



We shouldn't have a mailing list where thousands of recipients can receive
large emails :-)



I'm eager to know why! Thank you very much!



Force fields, including water models, are parameterized to reproduce
certain theoretical or experimental data. To the extent that those targets
resemble other observables, you might expect agreement. You should expect
the level of such agreement to differ between models.



Especially when the methods applied are different.  The original Berendsen 1987 
SPC/E paper used very different methods relative to what the OP shows in the 
.mdp file.  In the original paper, the cutoffs were 0.9 nm, no PME, Berendsen's 
weak coupling method for temperature and pressure, temperature at 300 K, SHAKE 
instead of SETTLE (doubt that will make much of a difference, honestly), and the 
simulations were only a whopping 27.5 ps.  So there's no doubt in my mind that 
one will get different results.  Also, every paper that involves SPC/E in some 
form of comparison reports a different diffusion coefficient (though this is not 
unique to SPC/E), and most are not corrected for system size effects. 
Reproducibility is a considerable challenge...


-Justin


Mark



TIP4P simulation results:
Statistics over 3001 steps using 31 frames

   Energies (kJ/mol)
LJ (SR)  Disper. corr.   Coulomb (SR)   Coul. recip.  Potential
  7.71821e+03   -2.05973e+02   -4.93956e+049.58136e+01
 -4.17875e+04
Kinetic En.   Total EnergyTemperature Pres. DC (bar) Pressure (bar)
7.30494e+03   -3.44826e+042.93006e+02   -1.14009e+021.41117e+00


SPC/E simulation results:

Statistics over 3001 steps using 31 frames


   Energies (kJ/mol)
LJ (SR)  Disper. corr.   Coulomb (SR)   Coul. recip.  Potential
9.23352e+03   -2.12206e+02   -5.62576e+049.63633e+01   -4.71400e+04
Kinetic En.   Total EnergyTemperature Pres. DC (bar) Pressure (bar)
7.30494e+03   -3.98350e+042.93005e+02   -1.18040e+021.64654e+00




title= OPLS MD waterethanol siimulation
; Run parameters
integrator= md; leap-frog integrator
nsteps= 3000; 1 * 2000 = 2 ps (20 ns)
dt= 0.001; 1 fs
; Output control
nstxout= 10; save coordinates every 10.0 ps
nstvout= 10; save velocities every 10.0 ps
nstenergy= 10; save energies every 10.0 ps
nstlog= 10; update log file every 10.0 ps
nstxout-compressed  = 10  ; save compressed coordinates every 10.0
ps
define  =-DEFLEXIBLE ; flexible
water
compressed-x-grps   = System; replaces xtc-grps
; Bond parameters
 continuation= yes; Restarting after NPT
 constraint_algorithm= lincs; holonomic constraints
 constraints= none; H bonds (even heavy atom-H bonds)
constrained
 lincs_iter= 1; accuracy of LINCS
 lincs_order= 4; also related to accuracy
; Neighborsearching
cutoff-scheme   = Verlet
ns_type= grid; search neighboring grid cells
nstlist= 10; 20 fs, largely irrelevant with Verlet scheme
rcoulomb= 1.2; short-range electrostatic cutoff (in nm)
rvdw= 1.2; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype= PME; Particle Mesh Ewald for long-range electrostatics
pme_order= 4; cubic interpolation
fourierspacing= 0.12; grid spacing for FFT
; Temperature coupling is on
tcoupl= V-rescale; modified Berendsen thermostat
tc-grps= SOL  ; two coupling groups - more accurate
tau_t= 0.1; time constant, in ps
ref_t= 293  ; reference temperature, one for each group, in K
; Pressure coupling is on
pcoupl= Parrinello-Rahman; Pressure coupling on in NPT
pcoupltype= isotropic; uniform scaling of box vectors
tau_p= 2.0; time constant, in ps
ref_p= 1.0; reference pressure, in bar
compressibility = 4.5e-5; isothermal compressibility of
water, bar^-1
; Periodic boundary conditions
pbc= xyz; 3-D PBC
; Dispersion 

Re: [gmx-users] assignment of atomname1 and atomname2 in param_example file in GridMAT-MD

[Justin Lemkul]

On 5/11/17 5:24 AM, manindersingh rajawat wrote:

Dear GROMACS users,

I have set up a system for membrane simulation. Membrane contains 103 POPC
and 25 POPS molecules. I am using GridMAT-MD for calculating area per lipid
and thickness of the membrane by following GridMAT-MD tutorial. But
confused how and what to assign in atomname1 and atomname2 in its
param_example file.



Use something that is representative of the lipid head group.  For most 
phospholipids, use the P atom.  If you have further questions, direct them 
off-list to me as this is not a GROMACS problem.


-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] Calculation of dielectric constant

[Saumyak Mukherjee]

Thanks a lot Philip for your suggestion.

On 11 May 2017 at 13:41, Philip Loche  wrote:

> Hey Saumyak,
>
> g_dipoles applies the formula derived by Neumann (
> http://dx.doi.org/10.1080/00268978300102721), which is proportional to
> 1/V (box volume) and can not applied in your case. According to your
> description you should consider to calculate a dielectric profile around
> your protein. If its spherical take a look at the following paper:
> http://link.aps.org/doi/10.1103/PhysRevE.92.032718 <
> http://link.aps.org/doi/10.1103/PhysRevE.92.032718>
>
> Best
>
> Philip
>
> —
>
> Philip Loche
> Department of Physics
> Freie Universität Berlin
> Arnimallee 14, Room 0.3.30
> 14195 Berlin, Germany
>
> phone: +49 (0)30/838-55279
>
> > On 11. May 2017, at 07:49, Saumyak Mukherjee <
> mukherjee.saumya...@gmail.com> wrote:
> >
> > Dear users,
> >
> > I am trying to calculate the dielectric constant of water within a
> certain
> > distance of protein.
> >
> > For this I have generated an ordered trajectory using trjorder. Then I
> have
> > created an index group for an estimated number of water molecules
> (average
> > of the numbers in nshell.xvg generated from trjorder).
> >
> > Thereafter, I have used the g_dipoles program along with the flag -eps to
> > get dielectric constant. But the results are incorrect, as I am expecting
> > the dielectric constant to be greater than 40 or so, but it is showing
> ~3.5.
> >
> > What am I doing wrong? Any help would be highly appreciated.
> >
> > Thanks and regards,
> > Saumyak
> > --
> > 
> > *Saumyak Mukherjee*
> >
> > Junior Research Fellow
> > Prof. Biman Bagchi's Group
> > Solid State and Structural Chemistry Unit
> > Indian Institute of Science
> > Bangalore - 560012
> >
> > Mob : 8017292426
> > Alternative e-mail : saumyakmukher...@gmail.com
> >smukher...@sscu.iisc.ernet.in
> > 
> > --
> > 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.




-- 

*Saumyak Mukherjee*

Junior Research Fellow
Prof. Biman Bagchi's Group
Solid State and Structural Chemistry Unit
Indian Institute of Science
Bangalore - 560012

Mob : 8017292426
Alternative e-mail : saumyakmukher...@gmail.com
smukher...@sscu.iisc.ernet.in

-- 
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] How to use "gmx view" on Ubuntu?

[ZHANG Cheng]

Dear Justin, Mark and Frank,
Thank you so much for explaining it. Now I think I understand it.


I was confusing because I thought, I should install/compile the xorg-dev by 
using "cmake". But actually, what you meant is, use "cmake" to compile Gromacs 
after installing xorg-dev.


I tried as Frank suggested for the cmake command line, including 
"-DGMX_X11=on", using my original Gromacs 5.0.4 files. But some problems 
occurred on the cmake. But I managed to use VMD to visualise my xtc file. So it 
is okay for me now.


Thank you for all your help!


Yours sincerely
Cheng




-- Original --
From:  "ZHANG Cheng";<272699...@qq.com>;
Date:  Thu, May 11, 2017 04:27 AM
To:  "gromacs.org_gmx-users"; 

Subject:  Re:   Re:   Re:  How to use "gmx view" on Ubuntu?



Dear Mark,
Yes, Justin gave me a link for Gromacs-2016.3, but not the link for xorg-dev.


Do you mean, I should uninstall my current Gromacs 5.0.4, then install the 
Gromacs-2016.3, so that I can use the "gmx view" in the Gromacs-2016.3 on my 
Ubuntu?


Cheng




-- Original --
From:  "ZHANG Cheng";<272699...@qq.com>;
Date:  Thu, May 11, 2017 02:50 AM
To:  "gromacs.org_gmx-users"; 

Subject:  Re:  Re:   Re:  How to use "gmx view" on Ubuntu?



Dear Justin,
Thank you for your link. I know how to install the Gromacs based on your link. 
But do you know where I can download the tar.gz file so that I can compile it 
and then use "gmx view"?


Thank you.


Cheng




-- Original --
From:  "ZHANG Cheng";<272699...@qq.com>;
Date:  Wed, May 10, 2017 11:27 PM
To:  "gromacs.org_gmx-users"; 

Subject:  Re:  Re:  How to use "gmx view" on Ubuntu?



Dear Mark Abraham,
Thank you very much for your updating. 


Sorry, could you please tell me where can I download the cmake file, or which 
command line to use, in order to compiling & installing? 


I have tried the below and they both worked, but "gmx view ??" still not work, 
with the same issue "Compiled without X-Windows - can not run viewer."
) sudo apt-get install xorg-dev
) sudo apt-get install xorg-dev libglu1-mesa-dev



(Sorry, I am not sure how to find the files to download)


Thank you.


Yours sincerely
Cheng


-- Original --
From:  "ZHANG Cheng";<272699...@qq.com>;
Date:  Wed, May 10, 2017 04:22 AM
To:  "gromacs.org_gmx-users"; 

Subject:  Re: How to use "gmx view" on Ubuntu?



Dear Gromacs and Mark Abraham,
Sorry, I am not familiar with compiling issue. 


I installed xorg-dev by:
sudo apt-get install xorg-dev


Then, I type 
cmake .. -DGMX_X11=on


But I was told: 
CMake Error: The source directory "/home/lanselibai/Cheng/gromacs/... ..." does 
not appear to contain CMakeLists.txt.


What I should do next?


Thank you.


Cheng






-- Original --
From:  "ZHANG Cheng";<272699...@qq.com>;
Date:  Wed, May 10, 2017 02:58 AM
To:  "gromacs.org_gmx-users"; 

Subject:  How to use "gmx view" on Ubuntu?



Dear Gromacs,
I am running Gromacs 5.0.4 on Ubuntu 14.04.1 LTS.


When I run "gmx view ..." as below:


gmx view -f dppc-md.xtc -s dppc-md.tpr


I got the following:


Compiled without X-Windows - can not run viewer.





Can I ask how to use "gmx view" on Ubuntu? Thank you.


Yours sincerely
Cheng
-- 
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] Using the md integrator for calculating free energy of solvation

[Dan Gil]

Thank you Dr. Lemkul,

The simulations ran, but I have some a question about the results.

Is it natural to have a nonzero mass_lambda dH/dL? My intuition was that
the potential energy does not depend on mass, and my thinking was that the
mass contribution should be very small, if not zero.

Best Regards,

Dan

On Wed, May 3, 2017 at 10:46 AM, Justin Lemkul  wrote:

>
>
> On 5/3/17 10:43 AM, Dan Gil wrote:
>
>> Hi Dr. Kausar,
>>
>> If I want the molecule to have all of its nonbonded interactions with the
>> rest of the system (solvent) at each lambda, don't I want to have vdwq
>> option ON for both couple-lambda0 and couple-lambda1?
>>
>> What I am attempting to do is mutate molecule A to molecule B through the
>> simulation.
>>
>>
> Relative free energy calculations are a little bit of a strange beast; you
> do need to set the B-state (couple-lambda1) to "none" and the B-state types
> and charges are read from [atoms].
>
>
> -Justin
>
> Best Regards,
>>
>> Dan
>>
>> On Wed, May 3, 2017 at 1:04 AM, Tasneem Kausar > >
>> wrote:
>>
>> In the last line of mdp file couple-lambda0 and couple-lambda1 have the
>>> same vdwq variable. It indicates A and B state of system. So grompp is
>>> complaining about the identical states.
>>>
>>> You can also see the alchemistry.org for detail information about the
>>> free
>>> energy calculations.
>>>
>>> On Wed, May 3, 2017 at 2:14 AM, Dan Gil  wrote:
>>>
>>> Thank you Dr. Lemkul,

 Following your advice, I was able to calculate the solvation free energy

>>> of
>>>
 two molecules that I am interested in. The results, unfortunately, were

>>> not
>>>
 what I expected. I want to try an alternative method, that is, the
 alchemical thermodynamic integration method.

 If I understand this correctly, I should be able to calculate the free
 energy difference between two systems A and B as I change the topology
 of
 the molecule from A to B in small steps.

 I wasn't able to find a tutorial online, but I attempted to try the

>>> method
>>>
 and I obtained this error from grompp:

 WARNING 1 [file grompp.mdp, line 43]:
   The lambda=0 and lambda=1 states for coupling are identical

 I was hoping if you could help me understand what I am missing in my
 topology or mdp file.

 Topology:
  [ moleculetype ]
 HEPT 3

  [ atoms ]
  1   opls_135  1   HEPTCH3  1 -0.180 12.011
  opls_9610.36012.011
  2   opls_136  1   HEPTCH2  2 -0.120 12.011
  opls_9620.24012.011
  3   opls_136  1   HEPTCH2  3 -0.120 12.011
  opls_9620.24012.011
  4   opls_136  1   HEPTCH2  4 -0.120 12.011
  opls_9620.24012.011
  5   opls_136  1   HEPTCH2  5 -0.120 12.011
  opls_9620.24012.011
  6   opls_136  1   HEPTCH2  6 -0.120 12.011
  opls_9620.24012.011
  7   opls_135  1   HEPTCH3  7 -0.180 12.011
  opls_9610.36012.011
  8   opls_140  1   HEPT  H  1  0.060  1.008
  opls_965   -0.12018.998
  9   opls_140  1   HEPT  H  1  0.060  1.008
  opls_965   -0.12018.998
 10   opls_140  1   HEPT  H  1  0.060  1.008
  opls_965   -0.12018.998
 11   opls_140  1   HEPT  H  2  0.060  1.008
  opls_965   -0.12018.998
 12   opls_140  1   HEPT  H  2  0.060  1.008
  opls_965   -0.12018.998
 13   opls_140  1   HEPT  H  3  0.060  1.008
  opls_965   -0.12018.998
 14   opls_140  1   HEPT  H  3  0.060  1.008
  opls_965   -0.12018.998
 15   opls_140  1   HEPT  H  4  0.060  1.008
  opls_965   -0.12018.998
 16   opls_140  1   HEPT  H  4  0.060  1.008
  opls_965   -0.12018.998
 17   opls_140  1   HEPT  H  5  0.060  1.008
  opls_965   -0.12018.998
 18   opls_140  1   HEPT  H  5  0.060  1.008
  opls_965   -0.12018.998
 19   opls_140  1   HEPT  H  6  0.060  1.008
  opls_965   -0.12018.998
 20   opls_140  1   HEPT  H  6  0.060  1.008
  opls_965   -0.12018.998
 21   opls_140  1   HEPT  H  7  0.060  1.008
  opls_965   -0.12018.998
 22   opls_140  1   HEPT  H  7  0.060  1.008
  opls_965   -0.12018.998
 23   opls_140  1   HEPT  H  7  0.060  1.008
 

Re: [gmx-users] ligand moving out during umbrella sampling

[abhisek Mondal]

On Thu, May 11, 2017 at 11:02 AM, abhisek Mondal 
wrote:

> Hi,
>
> Thank you for the explanation. It really cleared some concepts. But I'm
> still having my ligand moving in this step. I have modified the code as:
> ; Pull code
> pull= umbrella
> pull_ngroups= 1
> pull_group0 = Protein_chain_A
> pull_group1 = ACO
> pull_geometry   = direction ; simple distance increase
> pull_dim   = Y Y Y ; not to allow ligand move along other
> dir
> pull_rate1 = 0.0
> pull_k1   = 1000   ; kJ mol^-1 nm^-2
> pull_start   = yes   ; define initial COM distance > 0
> pull_vec1   = 0 0 -1
>
> The ligand was previously moving along x,y direction when I was using
> pull_dim  = N N Y. So I changed it to Y in all direction and provided 0 as
> vector  and pull_rate1=0.0, so that it does not move much. But at the end
> of a 10ns run, I see that the ligand is still moving as it was earlier.
>
It shows me:
 Pull group  natoms  pbc atom  distance at start reference at t=0
   0  1132936665
   159  1618  -1.555-1.555
Is it ok withe negative value ? Anyway this setup is not working.

>
> I have choose my reaction coordinate to be along -Z axis and want to apply
> biasing potential accordingly with restraining the ligand movement. Can you
> please suggest where am I failing with this code ?
>
> Thank you.
>
>
>
> On Tue, May 9, 2017 at 1:11 AM, Justin Lemkul  wrote:
>
>>
>>
>> On 5/8/17 10:00 AM, abhisek Mondal wrote:
>>
>>> On Sun, May 7, 2017 at 11:37 PM, Justin Lemkul  wrote:
>>>
>>>

 On 5/7/17 1:57 AM, abhisek Mondal wrote:

 Hi,
>
> For your ease of understanding regarding what is happening during this
> above said umbrella-mdrun, I have shared the trajectory video file the
> following link.
> https://drive.google.com/drive/folders/0B6O-L5Y7BiGJQ1FIc2tIRFE2dE0
>
> Is this normal given that the mdp code being used ? I basically have no
> idea with this step, so please help me out. I'm using gromacs-4.6.2.
>
>
> Your setup is incorrect.  You're applying a biasing potential only
 along
 z, so the ligand can move freely along x and y.  A protein-ligand
 complex
 has spherical symmetry, so you should set the reaction coordinate to the
 vector connecting the ligand with some suitable subset of interacting
 protein residues.

>>>
>>>
>>> I don't get it.
>>> We are trying not to move our configuration (generated after pulling
>>> simulation) along the reaction coordinate, so for restraining we are
>>> supposed to set pull_rate1=0.0.
>>>
>>
>> Of course.  But you said set pull_k = 0 which does not make sense.  The
>> pulling rate *is* zero during umbrella sampling (no net displacement,
>> restrain to the specified distance along the reaction coordinate) and
>> pulling force constant should be non-zero.
>>
>> If applying biasing potential only along z is causing movement along x and
>>> y then what if we apply the biasing potential along x,y,z ? Will it cause
>>> any good in restraining the ligand?
>>>
>>>
>> This is how it should be done.  The reaction coordinate should be
>> suitably defined based on the geometry of the system.  As I suggested
>> before, choose some representative residues in the active site as one group
>> and the ligand as the other.  Thus defines the reaction coordinate without
>> any presupposition of anything being aligned with a Cartesian axis, which
>> is rarely the case.
>>
>> Moreover, you said previously "A protein-ligand complex has spherical
>>> symmetry, so you should set the reaction coordinate to the vector
>>> connecting the ligand with some suitable subset of interacting protein
>>> residues.". It is really unclear to me, could you please give me some
>>> examples to understand it more simply? I had pulled the ligand along -z
>>> axis, doesn't it mean that the reaction coordinate is to be that way ?
>>> The
>>> fact that I'm struggling with is to restrain the pull configurations for
>>> further sampling.
>>>
>>>
>> The reaction coordinate is whatever you define it to be.  Whether or not
>> pulling along the z-axis makes sense depends on the orientation of the
>> system and the intrinsic geometry.  In your case, it doesn't make sense.
>> In my case (the tutorial, the unidirectional growth of an amyloid fibril)
>> it does make sense to use a single Cartesian axis for the SMD portion and
>> subsequent umbrella sampling.
>>
>> I'm really a beginner, so maybe I'm asking stupid questions. Please give
>>> me
>>> some advise. I'm really unable to decipher the scenario in comparison to
>>> your amyloid article in JPCB.
>>>
>>>
>> You should read the article to understand why I did what I did in the
>> tutorial, and then move on to reading articles that are more similar to
>> your case.  These will 

Re: [gmx-users] Running Gromacs on GPUs

[Mark Abraham]

Hi,

Each of the four simulations on your node has its own particle-particle
rank, and we knew that having multiple of them didn't work in 5.1.x, so
disabled it. But it was fixed for the 2016 release.

Mark

On Thu, 11 May 2017 13:00 Giannis Gl  wrote:

> Dear Gromacs users,
>
> I am trying to run a multiple walkers metadynamics simulation on Gromacs
> 5.1.4 using a machine that has 12 CPUs and 4 GPUs.
> I have compiled Gromacs using the following schems:
>
> cmake .. -DGMX_BUILD_OWN_FFTW=ON -DREGRESSIONTEST_DOWNLOAD=ON
> -DGMX_MPI=on -DGMX_GPU=on -DGMX_USE_OPENCL=on
>
> and the command that I use to run the simulation is a follows:
>
> mpirun -np 4 gmx_mpi mdrun -s topol -multi 4 -gpu_id 0123 -x traj
>
> However, I get the following error with or without the plumed flag:
>
> Running on 1 node with total 6 cores, 12 logical cores, 4 compatible GPUs
> ---
> Program gmx mdrun, VERSION 5.1.4
> Source code file:
> /usr/local/gromacs-5.1.4/src/gromacs/gmxlib/gmx_detect_hardware.cpp, line:
> 1203
>
> Fatal error:
> The OpenCL implementation only supports using exactly one PP rank per node
> For more information and tips for troubleshooting, please check the GROMACS
> website at http://www.gromacs.org/Documentation/Errors
> ---
>
> Halting parallel program gmx mdrun on rank 1 out of 4
> --
> MPI_ABORT was invoked on rank 1 in communicator MPI_COMM_WORLD
> with errorcode 1.
>
> NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
> You may or may not see output from other processes, depending on
> exactly when Open MPI kills them.
> --
>
> I tried the following command but still I get the same error,
>
> mpirun -np 4 gmx_mpi mdrun -s topol -multi 4 -gpu_id 0123 -ntomp 1
>
> Did anyone have the same issue at some point and if yes, how did you solve
> it?
>
> Best wishes,
> Yiannis
> --
> 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.

[gmx-users] Running Gromacs on GPUs

[Giannis Gl]

Dear Gromacs users,

I am trying to run a multiple walkers metadynamics simulation on Gromacs
5.1.4 using a machine that has 12 CPUs and 4 GPUs.
I have compiled Gromacs using the following schems:

cmake .. -DGMX_BUILD_OWN_FFTW=ON -DREGRESSIONTEST_DOWNLOAD=ON
-DGMX_MPI=on -DGMX_GPU=on -DGMX_USE_OPENCL=on

and the command that I use to run the simulation is a follows:

mpirun -np 4 gmx_mpi mdrun -s topol -multi 4 -gpu_id 0123 -x traj

However, I get the following error with or without the plumed flag:

Running on 1 node with total 6 cores, 12 logical cores, 4 compatible GPUs
---
Program gmx mdrun, VERSION 5.1.4
Source code file:
/usr/local/gromacs-5.1.4/src/gromacs/gmxlib/gmx_detect_hardware.cpp, line:
1203

Fatal error:
The OpenCL implementation only supports using exactly one PP rank per node
For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors
---

Halting parallel program gmx mdrun on rank 1 out of 4
--
MPI_ABORT was invoked on rank 1 in communicator MPI_COMM_WORLD
with errorcode 1.

NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.
--

I tried the following command but still I get the same error,

mpirun -np 4 gmx_mpi mdrun -s topol -multi 4 -gpu_id 0123 -ntomp 1

Did anyone have the same issue at some point and if yes, how did you solve
it?

Best wishes,
Yiannis
-- 
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] assignment of atomname1 and atomname2 in param_example file in GridMAT-MD

[manindersingh rajawat]

Dear GROMACS users,

I have set up a system for membrane simulation. Membrane contains 103 POPC
and 25 POPS molecules. I am using GridMAT-MD for calculating area per lipid
and thickness of the membrane by following GridMAT-MD tutorial. But
confused how and what to assign in atomname1 and atomname2 in its
param_example file.

Thanks and regards,
Maninder
-- 
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] Calculation of dielectric constant

[Philip Loche]

Hey Saumyak,

g_dipoles applies the formula derived by Neumann 
(http://dx.doi.org/10.1080/00268978300102721), which is proportional to 1/V 
(box volume) and can not applied in your case. According to your description 
you should consider to calculate a dielectric profile around your protein. If 
its spherical take a look at the following paper: 
http://link.aps.org/doi/10.1103/PhysRevE.92.032718 


Best

Philip

—

Philip Loche
Department of Physics
Freie Universität Berlin
Arnimallee 14, Room 0.3.30
14195 Berlin, Germany

phone: +49 (0)30/838-55279

> On 11. May 2017, at 07:49, Saumyak Mukherjee  
> wrote:
> 
> Dear users,
> 
> I am trying to calculate the dielectric constant of water within a certain
> distance of protein.
> 
> For this I have generated an ordered trajectory using trjorder. Then I have
> created an index group for an estimated number of water molecules (average
> of the numbers in nshell.xvg generated from trjorder).
> 
> Thereafter, I have used the g_dipoles program along with the flag -eps to
> get dielectric constant. But the results are incorrect, as I am expecting
> the dielectric constant to be greater than 40 or so, but it is showing ~3.5.
> 
> What am I doing wrong? Any help would be highly appreciated.
> 
> Thanks and regards,
> Saumyak
> -- 
> 
> *Saumyak Mukherjee*
> 
> Junior Research Fellow
> Prof. Biman Bagchi's Group
> Solid State and Structural Chemistry Unit
> Indian Institute of Science
> Bangalore - 560012
> 
> Mob : 8017292426
> Alternative e-mail : saumyakmukher...@gmail.com
>smukher...@sscu.iisc.ernet.in
> 
> -- 
> 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.

[gmx-users] Open GROMACS researcher position in Stockholm

[Erik Lindahl]

Hi,


Normally we don’t allow job announcements on these lists, with the one
exception for jobs that are directly related to GROMACS developement - and
I’m happy to announce that we have a new opening at the main GROMACS site
in Stockholm with application deadline June 16!


In particular, this is related to the new EU-funded BioExcel
Center-of-Excellence, and in this case we are looking for somebody who
would love to work with outreach, improving the way simulations work from a
user point-of-view (rather than only bottom-up coding), designing better
ways to solve concrete biophysics problems with GROMACS, and not least with
an interest in helping with arranging outreach, tutorials and documentation.


Basically, we are looking for somebody who loves all parts of simulations,
research, and computing - and who wants to join the team. There will be
plenty of opportunities to engage in cutting-edge research as part of the
work, and there are a number of international collaborations.


If you have any questions, contact either me or Mark Abraham <
mark.j.abra...@gmail.com>, and do help us with spreading the information
about the position:


https://www.kth.se/en/om/work-at-kth/lediga-jobb/what:job/jobID:150746/


All the best,


Erik

--
Erik Lindahl 
Professor of Biophysics, Dept. Biochemistry & Biophysics, Stockholm
University
Science for Life Laboratory, Box 1031, 17121 Solna, Sweden
-- 
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] soft "exit"

[Daskalakis Vangelis]

Thank you Mark.
That's very helpful
Vangelis


On Thu, May 11, 2017 at 10:00 AM, Mark Abraham 
wrote:

> Hi,
>
> On Thu, May 11, 2017 at 8:54 AM Daskalakis Vangelis 
> wrote:
>
> > Dear all,
> > I am running Gromacs 5.1.4 on a cluster with SLURM workload manager. How
> > can I terminate the running gromacs jobs (i.e. send the TERM/ INT
> signals)
> > but also allow time for each job to "soft exit" on the next step ? By
> "soft
> > exit" I mean that the job should be like it is finished (all the usual
> > output will be written to file), but before its end simulation time
> defined
> > in the tpr. There is of course the "scancel" command within slurm, but
> this
> > does not guarantee that the job will have the adequate time to "soft
> exit".
> >
>
> scancel -s SIGINT
>
> directs slurm to give the processes in the batch job that signal. When
> GROMACS gets it, it writes a checkpoint and does an orderly exit at the
> next neighbour-search step.
>
> Using mdrun -maxh 24 is a useful way to have mdrun know in advance that it
> should stop before 24 hours, which might be when a job scheduler starts
> sending signals.
>
> Mark
>
> Thanks,
> > Vangelis.
> > 
> >
> > 
> > --
> > 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.
>
-- 
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] soft "exit"

[Mark Abraham]

Hi,

On Thu, May 11, 2017 at 8:54 AM Daskalakis Vangelis 
wrote:

> Dear all,
> I am running Gromacs 5.1.4 on a cluster with SLURM workload manager. How
> can I terminate the running gromacs jobs (i.e. send the TERM/ INT signals)
> but also allow time for each job to "soft exit" on the next step ? By "soft
> exit" I mean that the job should be like it is finished (all the usual
> output will be written to file), but before its end simulation time defined
> in the tpr. There is of course the "scancel" command within slurm, but this
> does not guarantee that the job will have the adequate time to "soft exit".
>

scancel -s SIGINT

directs slurm to give the processes in the batch job that signal. When
GROMACS gets it, it writes a checkpoint and does an orderly exit at the
next neighbour-search step.

Using mdrun -maxh 24 is a useful way to have mdrun know in advance that it
should stop before 24 hours, which might be when a job scheduler starts
sending signals.

Mark

Thanks,
> Vangelis.
> 
>
> 
> --
> 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.

[gmx-users] soft "exit"

[Daskalakis Vangelis]

Dear all,
I am running Gromacs 5.1.4 on a cluster with SLURM workload manager. How
can I terminate the running gromacs jobs (i.e. send the TERM/ INT signals)
but also allow time for each job to "soft exit" on the next step ? By "soft
exit" I mean that the job should be like it is finished (all the usual
output will be written to file), but before its end simulation time defined
in the tpr. There is of course the "scancel" command within slurm, but this
does not guarantee that the job will have the adequate time to "soft exit".

Thanks,
Vangelis.



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