[gmx-users] Regarding generating topology file.

[Ravi Kumar Venkatraman]

Dear all,
 I have got some topology and parameter file for some small
organic molecules from swissparam in *.mol2 format. I am trying to generate
*.top, *.gro & *.itp files using topolbuild. But I dont know how to do. So
can anybody give me a link or some materials related to this.

Thank you.

*With Regards,**
*
*
*
*Ravi Kumar Venkatraman,*
*c/o Prof. Siva Umapathy,**
*
*IPC Dept., IISc.,*
*Bangalore, INDIA.*
*Phone No: +91-9686933963*
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Regarding generating topology file.

[Mark Abraham]

On 02/08/11, Ravi Kumar Venkatraman   wrote:
> Dear all,
>  I have got some topology and parameter file for some small 
> organic molecules from swissparam in *.mol2 format. I am trying to generate 
> *.top, *.gro & *.itp files using topolbuild. But I dont know how to do. So 
> can anybody give me a link or some materials related to this.
> 

Swissparam is supposed to provide these files, and have a tutorial for how to 
handle the subsequent details. Please search for these.

Mark
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] atom types

[Gavin Melaugh]

Hi Justin

Again thanks for the reply. I am not disagreeing with you but If I don't
include a [pairs] directive in the topology file (with gen_pairs =yes),
then there are no 1-4 LJ nor 1-4 Coulombic energies written in the log
file. When I include the [pair s] directive then both types of
interaction are written to the log file. Therefore does gen_pairs= yes +
[pairs] directive generate 1,4 LJ and 1,4 Coulomb according to fudge LJ
and QQ?

Thanks

Gavin

Justin A. Lemkul wrote:
>
>
> Gavin Melaugh wrote:
>> Hi Justin
>>
>> Sorry to labour on this but:
>>
>> I don't quite understand what you mean when you say that nonbonded pair
>> interactions are not Coulombic. Surely nonbonded charged atoms interact
>> with each other, when close enough? (or by Nonbonded pair interactions
>> do you explcitly mean 1-4,1-5, etc.)
>> What are the 1-4 Coulombic interactions generated by, if not by
>> gen_pairs =yes in my case?
>> I have read this section of the manual loads and though I had a
>> comprehensive understanding of it, but now I am confused again.
>>
>
> Sure, there are nonbonded interactions for 1-4, 1-5, etc.  But the
> purpose of [pairtype] generation is for LJ terms only.  They are
> special 1-4 interactions between different atomtypes.  Look at any
> force field for which [pairtypes] are listed - they have only C6 and
> C12 terms.  Charges are not used for these calculations, but they are
> applied later during the MD using normal Coulombic equations and FudgeQQ.
>
> -Justin
>
>> Thanks
>>
>> Gavin
>> Justin A. Lemkul wrote:
>>>
>>> Gavin Melaugh wrote:
 Hi Justin

 I have checked the tpr file. Now it seems to assign the the two
 type of
 CHs as the same atom type, but at the same time with the specified
 charge from the [atoms] directive, as I expected. Concerning 1-4
 interactions and gen_pairs =yes, my concern is this; from the pair
 list
 and using gen_pairs = yes, does grompp then take the 1-4 coulombic
 interaction for CH from the [atomtypes] directive (as the meaning of
 gen_pairs =yes)?
 Or does it assign the charge based on the atom index in the pair list?

>>> Charges are irrelevant for generation of pair interactions.  Nonbonded
>>> pair interactions are LJ, not Coulombic.  You will certainly have 1-4
>>> Coulombic interactions, but they are not generated by gen_pairs.  See
>>> manual section 5.3.4.
>>>
>>> -Justin
>>>
 Many Thanks

 Gavin

 Justin A. Lemkul wrote:
> Gavin Melaugh wrote:
>> Hi all
>>
>> A very quick question. I have an atom-type labelled CH in the
>> atom-types
>> with a particular charge, and in the atom list I assign some of
>> these
>> specific atoms with zero charge as below. When I generate 1,4
>> interactions using gen_pairs =yes, what charge for the CH type
>> does it
>> use? Does gromacs assign the CH with the different charge as a new
>> atom
>> type.
>>
> Charges set in [atomtypes] are not used.  The zero charge is
> assigned.  Verify this by using gmxdump on your .tpr file.
>
> -Justin
>
>> ;Parameter level
>> [defaults]
>> ; nbfunccomb-rule gen-pairsfudgeLJ fudgeQQ
>>  1 3  yes0.5 0.5
>>
>> [atomtypes]
>> ;type mass   charge  ptype sigma(nm) 
>> epsilon(kjmol-1)
>>CB 12.011000  0.00   A  0.355000 
>> 0.292880
>>CA 12.011000 -0.115000   A  0.355000 
>> 0.292880
>>HC  1.008000  0.115000   A  0.242000 
>> 0.125520
>>CU 13.019000  0.265000   A  0.35 
>> 0.334720
>>NU 14.007000 -0.597000   A  0.325000 
>> 0.711280
>>CH 13.019000  0.332000   A  0.385000 
>> 0.334720
>>C3 15.035000  0.00   A  0.391000 
>> 0.669440
>>C2 14.027000  0.00   A  0.390500 
>> 0.493712
>>
>> ;Molecular level
>> [moleculetype]
>> ;   name nrexcl
>> isotridecylcage  3
>>
>> [atoms]
>> .
>>  72  CH   1   CGECH  24  0.3320 13.0190
>>73  C2   1   CGEC2  25  0.
>> 14.0270
>>74  C2   1   CGEC2  25  0.
>> 14.0270
>>75  C2   1   CGEC2  25  0.
>> 14.0270
>>76  C2   1   CGEC2  26  0.
>> 14.0270
>>77  C2   1   CGEC2  26  0.
>> 14.0270
>>78  C2   1   CGEC2  26  0.
>> 14.0270
>>79  C2   1   CGEC2  27  0.
>> 14.0270
>>80  C2   1   CGEC2  27  0.
>>

Re: [gmx-users] atom types

[Mark Abraham]

On 02/08/11, Gavin Melaugh   wrote:
> Hi Justin
> 
> Again thanks for the reply. I am not disagreeing with you but If I don't
> include a [pairs] directive in the topology file (with gen_pairs =yes),
> then there are no 1-4 LJ nor 1-4 Coulombic energies written in the log
> file. When I include the [pair s] directive then both types of
> interaction are written to the log file. Therefore does gen_pairs= yes +
> [pairs] directive generate 1,4 LJ and 1,4 Coulomb according to fudge LJ
> and QQ?
> 

Does manual section 5.3.4 answer your question?

Mark


> 
> Thanks
> 
> Gavin
> 
> Justin A. Lemkul wrote:
> >
> >
> > Gavin Melaugh wrote:
> >> Hi Justin
> >>
> >> Sorry to labour on this but:
> >>
> >> I don't quite understand what you mean when you say that nonbonded pair
> >> interactions are not Coulombic. Surely nonbonded charged atoms interact
> >> with each other, when close enough? (or by Nonbonded pair interactions
> >> do you explcitly mean 1-4,1-5, etc.)
> >> What are the 1-4 Coulombic interactions generated by, if not by
> >> gen_pairs =yes in my case?
> >> I have read this section of the manual loads and though I had a
> >> comprehensive understanding of it, but now I am confused again.
> >>
> >
> > Sure, there are nonbonded interactions for 1-4, 1-5, etc.  But the
> > purpose of [pairtype] generation is for LJ terms only.  They are
> > special 1-4 interactions between different atomtypes.  Look at any
> > force field for which [pairtypes] are listed - they have only C6 and
> > C12 terms.  Charges are not used for these calculations, but they are
> > applied later during the MD using normal Coulombic equations and FudgeQQ.
> >
> > -Justin
> >
> >> Thanks
> >>
> >> Gavin
> >> Justin A. Lemkul wrote:
> >>>
> >>> Gavin Melaugh wrote:
>  Hi Justin
> 
>  I have checked the tpr file. Now it seems to assign the the two
>  type of
>  CHs as the same atom type, but at the same time with the specified
>  charge from the [atoms] directive, as I expected. Concerning 1-4
>  interactions and gen_pairs =yes, my concern is this; from the pair
>  list
>  and using gen_pairs = yes, does grompp then take the 1-4 coulombic
>  interaction for CH from the [atomtypes] directive (as the meaning of
>  gen_pairs =yes)?
>  Or does it assign the charge based on the atom index in the pair list?
> 
> >>> Charges are irrelevant for generation of pair interactions.  Nonbonded
> >>> pair interactions are LJ, not Coulombic.  You will certainly have 1-4
> >>> Coulombic interactions, but they are not generated by gen_pairs.  See
> >>> manual section 5.3.4.
> >>>
> >>> -Justin
> >>>
>  Many Thanks
> 
>  Gavin
> 
>  Justin A. Lemkul wrote:
> > Gavin Melaugh wrote:
> >> Hi all
> >>
> >> A very quick question. I have an atom-type labelled CH in the
> >> atom-types
> >> with a particular charge, and in the atom list I assign some of
> >> these
> >> specific atoms with zero charge as below. When I generate 1,4
> >> interactions using gen_pairs =yes, what charge for the CH type
> >> does it
> >> use? Does gromacs assign the CH with the different charge as a new
> >> atom
> >> type.
> >>
> > Charges set in [atomtypes] are not used.  The zero charge is
> > assigned.  Verify this by using gmxdump on your .tpr file.
> >
> > -Justin
> >
> >> ;Parameter level
> >> [defaults]
> >> ; nbfunc    comb-rule gen-pairs    fudgeLJ fudgeQQ
> >>  1 3  yes    0.5 0.5
> >>
> >> [atomtypes]
> >> ;type mass   charge  ptype sigma(nm) 
> >> epsilon(kjmol-1)
> >>    CB 12.011000  0.00   A  0.355000 
> >> 0.292880
> >>    CA 12.011000 -0.115000   A  0.355000 
> >> 0.292880
> >>    HC  1.008000  0.115000   A  0.242000 
> >> 0.125520
> >>    CU 13.019000  0.265000   A  0.35 
> >> 0.334720
> >>    NU 14.007000 -0.597000   A  0.325000 
> >> 0.711280
> >>    CH 13.019000  0.332000   A  0.385000 
> >> 0.334720
> >>    C3 15.035000  0.00   A  0.391000 
> >> 0.669440
> >>    C2 14.027000  0.00   A  0.390500 
> >> 0.493712
> >>
> >> ;Molecular level
> >> [moleculetype]
> >> ;   name nrexcl
> >> isotridecylcage  3
> >>
> >> [atoms]
> >> .
> >>  72  CH   1   CGE    CH  24  0.3320 13.0190
> >>    73  C2   1   CGE    C2  25  0.    
> >> 14.0270
> >>    74  C2   1   CGE    C2  25  0.    
> >> 14.0270
> >>    75  C2   1   CGE    C2  25  0.    
> >> 14.0270
> >>    76  C2   1   CGE    C2   

Re: [gmx-users] atom types

[Gavin Melaugh]

Hi Mark


Thanks for the reply.
I am currently reading that section of the manual and, unless I am
completely mistaken, it seems to vindicate what I am saying.
"Extra Lennard-Jones and electrostatic interactions between pairs of
atoms in a molecule can be added in the [pairs] section of a molecule
definition".
In my [atom types] directive I have atomtype, charge mass, sigma and
epsilon etc. All nonbonding parameters are then calculated according to
the combination rule (in my case 3). 1-4 interactions are then
calculated based on the information in [pairs] directive (all atoms are
three bond away). I just have the atom indices of each pair in this
directive therefore with gen_pairs = yes, the interaction parameters
between each pair (which are 1-4) are calculated based on Fudge LJ and
Fudge QQ (which are both 0.5 in my case). All of this in conjunction
with nrexcl =3.  Or am I completely wrong?
In my set up then, are 1-4 Coulombic interactions determined by the pair
list and fudge QQ?

Many Thanks

Gavin


Mark Abraham wrote:
>
>
> On 02/08/11, *Gavin Melaugh *  wrote:
>> Hi Justin
>>
>> Again thanks for the reply. I am not disagreeing with you but If I don't
>> include a [pairs] directive in the topology file (with gen_pairs =yes),
>> then there are no 1-4 LJ nor 1-4 Coulombic energies written in the log
>> file. When I include the [pair s] directive then both types of
>> interaction are written to the log file. Therefore does gen_pairs= yes +
>> [pairs] directive generate 1,4 LJ and 1,4 Coulomb according to fudge LJ
>> and QQ?
>
> Does manual section 5.3.4 answer your question?
>
> Mark
>
>>
>> Thanks
>>
>> Gavin
>>
>> Justin A. Lemkul wrote:
>> >
>> >
>> > Gavin Melaugh wrote:
>> >> Hi Justin
>> >>
>> >> Sorry to labour on this but:
>> >>
>> >> I don't quite understand what you mean when you say that nonbonded
>> pair
>> >> interactions are not Coulombic. Surely nonbonded charged atoms
>> interact
>> >> with each other, when close enough? (or by Nonbonded pair interactions
>> >> do you explcitly mean 1-4,1-5, etc.)
>> >> What are the 1-4 Coulombic interactions generated by, if not by
>> >> gen_pairs =yes in my case?
>> >> I have read this section of the manual loads and though I had a
>> >> comprehensive understanding of it, but now I am confused again.
>> >>
>> >
>> > Sure, there are nonbonded interactions for 1-4, 1-5, etc.  But the
>> > purpose of [pairtype] generation is for LJ terms only.  They are
>> > special 1-4 interactions between different atomtypes.  Look at any
>> > force field for which [pairtypes] are listed - they have only C6 and
>> > C12 terms.  Charges are not used for these calculations, but they are
>> > applied later during the MD using normal Coulombic equations and
>> FudgeQQ.
>> >
>> > -Justin
>> >
>> >> Thanks
>> >>
>> >> Gavin
>> >> Justin A. Lemkul wrote:
>> >>>
>> >>> Gavin Melaugh wrote:
>>  Hi Justin
>> 
>>  I have checked the tpr file. Now it seems to assign the the two
>>  type of
>>  CHs as the same atom type, but at the same time with the specified
>>  charge from the [atoms] directive, as I expected. Concerning 1-4
>>  interactions and gen_pairs =yes, my concern is this; from the pair
>>  list
>>  and using gen_pairs = yes, does grompp then take the 1-4 coulombic
>>  interaction for CH from the [atomtypes] directive (as the meaning of
>>  gen_pairs =yes)?
>>  Or does it assign the charge based on the atom index in the pair
>> list?
>> 
>> >>> Charges are irrelevant for generation of pair interactions. 
>> Nonbonded
>> >>> pair interactions are LJ, not Coulombic.  You will certainly have 1-4
>> >>> Coulombic interactions, but they are not generated by gen_pairs.  See
>> >>> manual section 5.3.4.
>> >>>
>> >>> -Justin
>> >>>
>>  Many Thanks
>> 
>>  Gavin
>> 
>>  Justin A. Lemkul wrote:
>> > Gavin Melaugh wrote:
>> >> Hi all
>> >>
>> >> A very quick question. I have an atom-type labelled CH in the
>> >> atom-types
>> >> with a particular charge, and in the atom list I assign some of
>> >> these
>> >> specific atoms with zero charge as below. When I generate 1,4
>> >> interactions using gen_pairs =yes, what charge for the CH type
>> >> does it
>> >> use? Does gromacs assign the CH with the different charge as a new
>> >> atom
>> >> type.
>> >>
>> > Charges set in [atomtypes] are not used.  The zero charge is
>> > assigned.  Verify this by using gmxdump on your .tpr file.
>> >
>> > -Justin
>> >
>> >> ;Parameter level
>> >> [defaults]
>> >> ; nbfunccomb-rule gen-pairsfudgeLJ fudgeQQ
>> >>  1 3  yes0.5 0.5
>> >>
>> >> [atomtypes]
>> >> ;type mass   charge  ptype sigma(nm)
>> >> epsilon(kjmol-1)
>> >>CB 12.011000  0.00   A  0.355000
>> >> 0.292880
>> >>CA 12.011000 

[gmx-users] Re: a question about order parameter of POPC

[Justin A. Lemkul]

There is a bug in g_order.  It produces -Scd values that are similar to 
experimental results and those from other simulations, but they are not correct, 
at least in the case of unsaturations.  Saturated carbons should be correct, if 
I recall.  So no, there is no satisfactory answer until the bug is fixed.  A 
code patch was proposed proposed some time ago, but it has not been officially 
implemented.


Please post anything further to gmx-users.

-Justin

Amir Mohsen Pourmousa wrote:

Hi Justin,

I just noticed on the mailing list of GROMACS that you had the same 
problem that I have now. The response that you got from one of those 
people was not satisfactory (I think it didn't satisfy you either):

http://www.mail-archive.com/gmx-users@gromacs.org/msg28555.html

Anyway I am trying to compute the order parameter for POPC. I followed 
the instruction on Gromacs manual regarding the double bonds but my 
graph is different from those in literature. Did you finally succeed in 
obtaining the correct order parameter for the oleoyl acyl chain of POPC 
which contains double bonds? If yes, would you tell me how?


Best regards,
Amir-Mohsen Pourmousa
 
---

Amir-Mohsen Pourmousa
PhD Candidate of Theoretical Physics,
Dept. of Applied Mathematics, University of Western Ontario
1151 Richmond St. North, London (ON),  Canada N6A 5B7
Tel: +1-519-661 2111 extension 88678
Cell: +1-519-854 3331
Email: apour...@uwo.ca 


--


Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.

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

Re: [gmx-users] atom types

[Mark Abraham]

On 02/08/11, Gavin Melaugh   wrote:
> Hi Mark
> 
> 
> Thanks for the reply.
> I am currently reading that section of the manual and, unless I am
> completely mistaken, it seems to vindicate what I am saying.
> "Extra Lennard-Jones and electrostatic interactions between pairs of
> atoms in a molecule can be added in the [pairs] section of a molecule
> definition".
> In my [atom types] directive I have atomtype, charge mass, sigma and
> epsilon etc. All nonbonding parameters are then calculated according to
> the combination rule (in my case 3). 1-4 interactions are then
> calculated based on the information in [pairs] directive (all atoms are
> three bond away). I just have the atom indices of each pair in this
> directive therefore with gen_pairs = yes, the interaction parameters
> between each pair (which are 1-4) are calculated based on Fudge LJ and
> Fudge QQ (which are both 0.5 in my case). All of this in conjunction
> with nrexcl =3.
> 

That will generate parameters for the interactions listed in [pairs] that do 
not have corresponding [pairtypes]. FudgeLJ and [nonbond_params] are used in 
such generation, per other parts of 5.7.


>   Or am I completely wrong?
> In my set up then, are 1-4 Coulombic interactions determined by the pair
> list and fudge QQ?
> 

If the contradiction you think exists is this one...

> > On 02/08/11, *Gavin Melaugh *  wrote:
> >> Hi Justin
> >>
> >> Again thanks for the reply. I am not disagreeing with you but If I don't
> >> include a [pairs] directive in the topology file (with gen_pairs =yes),
> >> then there are no 1-4 LJ nor 1-4 Coulombic energies written in the log
> >> file. When I include the [pair s] directive then both types of
> >> interaction are written to the log file. Therefore does gen_pairs= yes +
> >> [pairs] directive generate 1,4 LJ and 1,4 Coulomb according to fudge LJ
> >> and QQ?
> 

... then 5.3.4 indicates that the presence of a [pairs] directive will generate 
the 1,4 output fields. The parameters for that output are taken from 
[pairtypes]. If gen-pairs=yes then the parameters are generated, else some 
warning/error occurs. The example in 5.7.1 has some more explanation about the 
use of the fudge parameters.

Mark
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] atom types

[Gavin Melaugh]

Yes I think the example vindicates what I am saying as well. I suppose I
the "contradiction" ( I'll call it the point of confusion) you refer to
is perhaps when Justin (who is always more than helpful) said that
"Charges are irrelevant for generation of pair interactions.  Nonbonded
pair interactions are LJ, not Coulombic.  You will certainly have 1-4
Coulombic interactions, but they are not generated by gen_pairs.  See
manual section 5.3.4."

My sequence of 1-4 interaction generation should go like this I suppose:

e.g.
[pairs]
3  6
no parameters present  therefore get from [pairtypes] directive.
no [pairtypes] directive therefore get from [non_bonded parameters]
directive as gen pairs = yes
again no [non_bonded parameters] directive.
Therefore generate 1,4 interaction parameters based on the normal sigma
and epsilon values (comb rule 3) present in [atomtypes] directive, in
accordance with fudge LJ and QQ.

My point is, then in conclusion, that in this way surely the 1,4
electrostatic interactions are determined by the pair list and in my
case gen_pairs = yes no?

Many Thanks

Gavin


Mark Abraham wrote:
> On 02/08/11, *Gavin Melaugh *  wrote:
>> Hi Mark
>>
>>
>> Thanks for the reply.
>> I am currently reading that section of the manual and, unless I am
>> completely mistaken, it seems to vindicate what I am saying.
>> "Extra Lennard-Jones and electrostatic interactions between pairs of
>> atoms in a molecule can be added in the [pairs] section of a molecule
>> definition".
>> In my [atom types] directive I have atomtype, charge mass, sigma and
>> epsilon etc. All nonbonding parameters are then calculated according to
>> the combination rule (in my case 3). 1-4 interactions are then
>> calculated based on the information in [pairs] directive (all atoms are
>> three bond away). I just have the atom indices of each pair in this
>> directive therefore with gen_pairs = yes, the interaction parameters
>> between each pair (which are 1-4) are calculated based on Fudge LJ and
>> Fudge QQ (which are both 0.5 in my case). All of this in conjunction
>> with nrexcl =3.
>
> That will generate parameters for the interactions listed in [pairs]
> that do not have corresponding [pairtypes]. FudgeLJ and
> [nonbond_params] are used in such generation, per other parts of 5.7.
>
>>   Or am I completely wrong?
>> In my set up then, are 1-4 Coulombic interactions determined by the pair
>> list and fudge QQ?
>
> If the contradiction you think exists is this one...
>> > On 02/08/11, *Gavin Melaugh *  wrote:
>> >> Hi Justin
>> >>
>> >> Again thanks for the reply. I am not disagreeing with you but If I
>> don't
>> >> include a [pairs] directive in the topology file (with gen_pairs
>> =yes),
>> >> then there are no 1-4 LJ nor 1-4 Coulombic energies written in the log
>> >> file. When I include the [pair s] directive then both types of
>> >> interaction are written to the log file. Therefore does gen_pairs=
>> yes +
>> >> [pairs] directive generate 1,4 LJ and 1,4 Coulomb according to
>> fudge LJ
>> >> and QQ?
>
> ... then 5.3.4 indicates that the presence of a [pairs] directive will
> generate the 1,4 output fields. The parameters for that output are
> taken from [pairtypes]. If gen-pairs=yes then the parameters are
> generated, else some warning/error occurs. The example in 5.7.1 has
> some more explanation about the use of the fudge parameters.
>
> Mark 

-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] atom types

[Justin A. Lemkul]

Gavin Melaugh wrote:

Yes I think the example vindicates what I am saying as well. I suppose I
the "contradiction" ( I'll call it the point of confusion) you refer to
is perhaps when Justin (who is always more than helpful) said that
"Charges are irrelevant for generation of pair interactions.  Nonbonded
pair interactions are LJ, not Coulombic.  You will certainly have 1-4
Coulombic interactions, but they are not generated by gen_pairs.  See
manual section 5.3.4."



Charges *are* irrelevant - the information is not used when generating pairs, 
which I thought was the original question.  The charge information is used 
during MD, when the pair list tells mdrun which atoms interact in what way. 
Then you get Coul. 1-4 terms.  Perhaps I missed your point, but this whole 
thread started as "which charge is used to generate pairs?"  The answer is still 
none.



My sequence of 1-4 interaction generation should go like this I suppose:

e.g.
[pairs]
3  6
no parameters present  therefore get from [pairtypes] directive.
no [pairtypes] directive therefore get from [non_bonded parameters]
directive as gen pairs = yes
again no [non_bonded parameters] directive.
Therefore generate 1,4 interaction parameters based on the normal sigma
and epsilon values (comb rule 3) present in [atomtypes] directive, in
accordance with fudge LJ and QQ.

My point is, then in conclusion, that in this way surely the 1,4
electrostatic interactions are determined by the pair list and in my
case gen_pairs = yes no?



Yes.

-Justin


Many Thanks

Gavin


Mark Abraham wrote:

On 02/08/11, *Gavin Melaugh *  wrote:

Hi Mark


Thanks for the reply.
I am currently reading that section of the manual and, unless I am
completely mistaken, it seems to vindicate what I am saying.
"Extra Lennard-Jones and electrostatic interactions between pairs of
atoms in a molecule can be added in the [pairs] section of a molecule
definition".
In my [atom types] directive I have atomtype, charge mass, sigma and
epsilon etc. All nonbonding parameters are then calculated according to
the combination rule (in my case 3). 1-4 interactions are then
calculated based on the information in [pairs] directive (all atoms are
three bond away). I just have the atom indices of each pair in this
directive therefore with gen_pairs = yes, the interaction parameters
between each pair (which are 1-4) are calculated based on Fudge LJ and
Fudge QQ (which are both 0.5 in my case). All of this in conjunction
with nrexcl =3.

That will generate parameters for the interactions listed in [pairs]
that do not have corresponding [pairtypes]. FudgeLJ and
[nonbond_params] are used in such generation, per other parts of 5.7.


  Or am I completely wrong?
In my set up then, are 1-4 Coulombic interactions determined by the pair
list and fudge QQ?

If the contradiction you think exists is this one...

On 02/08/11, *Gavin Melaugh *  wrote:

Hi Justin

Again thanks for the reply. I am not disagreeing with you but If I

don't

include a [pairs] directive in the topology file (with gen_pairs

=yes),

then there are no 1-4 LJ nor 1-4 Coulombic energies written in the log
file. When I include the [pair s] directive then both types of
interaction are written to the log file. Therefore does gen_pairs=

yes +

[pairs] directive generate 1,4 LJ and 1,4 Coulomb according to

fudge LJ

and QQ?

... then 5.3.4 indicates that the presence of a [pairs] directive will
generate the 1,4 output fields. The parameters for that output are
taken from [pairtypes]. If gen-pairs=yes then the parameters are
generated, else some warning/error occurs. The example in 5.7.1 has
some more explanation about the use of the fudge parameters.

Mark 




--


Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.

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

[gmx-users] to clear up the confusion regarding fudge LJ and QQ

[Gavin Melaugh]

Hi Justin and Mark

Fudge LJ is only used when gen_pairs =yes. Fudge QQ is used regardless.
I suppose a more concise and sensible question is the following:
Are the atoms involved in the 1,4 Coulombic interactions specified in
the [pairs] directive? If not how else are they generated.

Cheers

Gavin

Gavin Melaugh wrote:
> Yes I think the example vindicates what I am saying as well. I suppose I
> the "contradiction" ( I'll call it the point of confusion) you refer to
> is perhaps when Justin (who is always more than helpful) said that
> "Charges are irrelevant for generation of pair interactions.  Nonbonded
> pair interactions are LJ, not Coulombic.  You will certainly have 1-4
> Coulombic interactions, but they are not generated by gen_pairs.  See
> manual section 5.3.4."
>
> My sequence of 1-4 interaction generation should go like this I suppose:
>
> e.g.
> [pairs]
> 3  6
> no parameters present  therefore get from [pairtypes] directive.
> no [pairtypes] directive therefore get from [non_bonded parameters]
> directive as gen pairs = yes
> again no [non_bonded parameters] directive.
> Therefore generate 1,4 interaction parameters based on the normal sigma
> and epsilon values (comb rule 3) present in [atomtypes] directive, in
> accordance with fudge LJ and QQ.
>
> My point is, then in conclusion, that in this way surely the 1,4
> electrostatic interactions are determined by the pair list and in my
> case gen_pairs = yes no?
>
> Many Thanks
>
> Gavin
>
>
> Mark Abraham wrote:
>   
>> On 02/08/11, *Gavin Melaugh *  wrote:
>> 
>>> Hi Mark
>>>
>>>
>>> Thanks for the reply.
>>> I am currently reading that section of the manual and, unless I am
>>> completely mistaken, it seems to vindicate what I am saying.
>>> "Extra Lennard-Jones and electrostatic interactions between pairs of
>>> atoms in a molecule can be added in the [pairs] section of a molecule
>>> definition".
>>> In my [atom types] directive I have atomtype, charge mass, sigma and
>>> epsilon etc. All nonbonding parameters are then calculated according to
>>> the combination rule (in my case 3). 1-4 interactions are then
>>> calculated based on the information in [pairs] directive (all atoms are
>>> three bond away). I just have the atom indices of each pair in this
>>> directive therefore with gen_pairs = yes, the interaction parameters
>>> between each pair (which are 1-4) are calculated based on Fudge LJ and
>>> Fudge QQ (which are both 0.5 in my case). All of this in conjunction
>>> with nrexcl =3.
>>>   
>> That will generate parameters for the interactions listed in [pairs]
>> that do not have corresponding [pairtypes]. FudgeLJ and
>> [nonbond_params] are used in such generation, per other parts of 5.7.
>>
>> 
>>>   Or am I completely wrong?
>>> In my set up then, are 1-4 Coulombic interactions determined by the pair
>>> list and fudge QQ?
>>>   
>> If the contradiction you think exists is this one...
>> 
 On 02/08/11, *Gavin Melaugh *  wrote:
 
> Hi Justin
>
> Again thanks for the reply. I am not disagreeing with you but If I
>   
>>> don't
>>>   
> include a [pairs] directive in the topology file (with gen_pairs
>   
>>> =yes),
>>>   
> then there are no 1-4 LJ nor 1-4 Coulombic energies written in the log
> file. When I include the [pair s] directive then both types of
> interaction are written to the log file. Therefore does gen_pairs=
>   
>>> yes +
>>>   
> [pairs] directive generate 1,4 LJ and 1,4 Coulomb according to
>   
>>> fudge LJ
>>>   
> and QQ?
>   
>> ... then 5.3.4 indicates that the presence of a [pairs] directive will
>> generate the 1,4 output fields. The parameters for that output are
>> taken from [pairtypes]. If gen-pairs=yes then the parameters are
>> generated, else some warning/error occurs. The example in 5.7.1 has
>> some more explanation about the use of the fudge parameters.
>>
>> Mark 
>> 
>
>   

-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] atom types

[Gavin Melaugh]

Hi Justin

I see that we may have got our wires crossed from the off.
Consider the [pairs] directive, which determines which atoms interact 
in a 1,4 manner.  Consider two atoms listed the [pairs] directive. From
the point of the Coulombic interaction between these two atoms I suppose
my original question should have been: Does mdrun, when calculating the
Coulombic potential between these two atoms, use the charges assigned to
the atoms in the [atomtypes] directive or [atoms] directive ?

Cheers

Gavin
 
Justin A. Lemkul wrote:
>
>
> Gavin Melaugh wrote:
>> Yes I think the example vindicates what I am saying as well. I suppose I
>> the "contradiction" ( I'll call it the point of confusion) you refer to
>> is perhaps when Justin (who is always more than helpful) said that
>> "Charges are irrelevant for generation of pair interactions.  Nonbonded
>> pair interactions are LJ, not Coulombic.  You will certainly have 1-4
>> Coulombic interactions, but they are not generated by gen_pairs.  See
>> manual section 5.3.4."
>>
>
> Charges *are* irrelevant - the information is not used when generating
> pairs, which I thought was the original question.  The charge
> information is used during MD, when the pair list tells mdrun which
> atoms interact in what way. Then you get Coul. 1-4 terms.  Perhaps I
> missed your point, but this whole thread started as "which charge is
> used to generate pairs?"  The answer is still none.
>
>> My sequence of 1-4 interaction generation should go like this I suppose:
>>
>> e.g.
>> [pairs]
>> 3  6
>> no parameters present  therefore get from [pairtypes] directive.
>> no [pairtypes] directive therefore get from [non_bonded parameters]
>> directive as gen pairs = yes
>> again no [non_bonded parameters] directive.
>> Therefore generate 1,4 interaction parameters based on the normal sigma
>> and epsilon values (comb rule 3) present in [atomtypes] directive, in
>> accordance with fudge LJ and QQ.
>>
>> My point is, then in conclusion, that in this way surely the 1,4
>> electrostatic interactions are determined by the pair list and in my
>> case gen_pairs = yes no?
>>
>
> Yes.
>
> -Justin
>
>> Many Thanks
>>
>> Gavin
>>
>>
>> Mark Abraham wrote:
>>> On 02/08/11, *Gavin Melaugh *  wrote:
 Hi Mark


 Thanks for the reply.
 I am currently reading that section of the manual and, unless I am
 completely mistaken, it seems to vindicate what I am saying.
 "Extra Lennard-Jones and electrostatic interactions between pairs of
 atoms in a molecule can be added in the [pairs] section of a molecule
 definition".
 In my [atom types] directive I have atomtype, charge mass, sigma and
 epsilon etc. All nonbonding parameters are then calculated
 according to
 the combination rule (in my case 3). 1-4 interactions are then
 calculated based on the information in [pairs] directive (all atoms
 are
 three bond away). I just have the atom indices of each pair in this
 directive therefore with gen_pairs = yes, the interaction parameters
 between each pair (which are 1-4) are calculated based on Fudge LJ and
 Fudge QQ (which are both 0.5 in my case). All of this in conjunction
 with nrexcl =3.
>>> That will generate parameters for the interactions listed in [pairs]
>>> that do not have corresponding [pairtypes]. FudgeLJ and
>>> [nonbond_params] are used in such generation, per other parts of 5.7.
>>>
   Or am I completely wrong?
 In my set up then, are 1-4 Coulombic interactions determined by the
 pair
 list and fudge QQ?
>>> If the contradiction you think exists is this one...
> On 02/08/11, *Gavin Melaugh *  wrote:
>> Hi Justin
>>
>> Again thanks for the reply. I am not disagreeing with you but If I
 don't
>> include a [pairs] directive in the topology file (with gen_pairs
 =yes),
>> then there are no 1-4 LJ nor 1-4 Coulombic energies written in
>> the log
>> file. When I include the [pair s] directive then both types of
>> interaction are written to the log file. Therefore does gen_pairs=
 yes +
>> [pairs] directive generate 1,4 LJ and 1,4 Coulomb according to
 fudge LJ
>> and QQ?
>>> ... then 5.3.4 indicates that the presence of a [pairs] directive will
>>> generate the 1,4 output fields. The parameters for that output are
>>> taken from [pairtypes]. If gen-pairs=yes then the parameters are
>>> generated, else some warning/error occurs. The example in 5.7.1 has
>>> some more explanation about the use of the fudge parameters.
>>>
>>> Mark 
>>
>

-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] atom types

[Mark Abraham]

On 02/08/11, Gavin Melaugh   wrote:
> Hi Justin
> 
> I see that we may have got our wires crossed from the off.
> Consider the [pairs] directive, which determines which atoms interact 
> in a 1,4 manner.  Consider two atoms listed the [pairs] directive. From
> the point of the Coulombic interaction between these two atoms I suppose
> my original question should have been: Does mdrun, when calculating the
> Coulombic potential between these two atoms, use the charges assigned to
> the atoms in the [atomtypes] directive or [atoms] directive ?
> 

The charges from [atoms] are used. The charges in [atomtypes] are never used in 
any forcefield currently used with GROMACS. It might not be possible for grompp 
to ever use them, but I'd have to check the code for that.

Mark


> 
> 
> Cheers
> 
> Gavin
>  
> Justin A. Lemkul wrote:
> >
> >
> > Gavin Melaugh wrote:
> >> Yes I think the example vindicates what I am saying as well. I suppose I
> >> the "contradiction" ( I'll call it the point of confusion) you refer to
> >> is perhaps when Justin (who is always more than helpful) said that
> >> "Charges are irrelevant for generation of pair interactions.  Nonbonded
> >> pair interactions are LJ, not Coulombic.  You will certainly have 1-4
> >> Coulombic interactions, but they are not generated by gen_pairs.  See
> >> manual section 5.3.4."
> >>
> >
> > Charges *are* irrelevant - the information is not used when generating
> > pairs, which I thought was the original question.  The charge
> > information is used during MD, when the pair list tells mdrun which
> > atoms interact in what way. Then you get Coul. 1-4 terms.  Perhaps I
> > missed your point, but this whole thread started as "which charge is
> > used to generate pairs?"  The answer is still none.
> >
> >> My sequence of 1-4 interaction generation should go like this I suppose:
> >>
> >> e.g.
> >> [pairs]
> >> 3  6
> >> no parameters present  therefore get from [pairtypes] directive.
> >> no [pairtypes] directive therefore get from [non_bonded parameters]
> >> directive as gen pairs = yes
> >> again no [non_bonded parameters] directive.
> >> Therefore generate 1,4 interaction parameters based on the normal sigma
> >> and epsilon values (comb rule 3) present in [atomtypes] directive, in
> >> accordance with fudge LJ and QQ.
> >>
> >> My point is, then in conclusion, that in this way surely the 1,4
> >> electrostatic interactions are determined by the pair list and in my
> >> case gen_pairs = yes no?
> >>
> >
> > Yes.
> >
> > -Justin
> >
> >> Many Thanks
> >>
> >> Gavin
> >>
> >>
> >> Mark Abraham wrote:
> >>> On 02/08/11, *Gavin Melaugh *  wrote:
>  Hi Mark
> 
> 
>  Thanks for the reply.
>  I am currently reading that section of the manual and, unless I am
>  completely mistaken, it seems to vindicate what I am saying.
>  "Extra Lennard-Jones and electrostatic interactions between pairs of
>  atoms in a molecule can be added in the [pairs] section of a molecule
>  definition".
>  In my [atom types] directive I have atomtype, charge mass, sigma and
>  epsilon etc. All nonbonding parameters are then calculated
>  according to
>  the combination rule (in my case 3). 1-4 interactions are then
>  calculated based on the information in [pairs] directive (all atoms
>  are
>  three bond away). I just have the atom indices of each pair in this
>  directive therefore with gen_pairs = yes, the interaction parameters
>  between each pair (which are 1-4) are calculated based on Fudge LJ and
>  Fudge QQ (which are both 0.5 in my case). All of this in conjunction
>  with nrexcl =3.
> >>> That will generate parameters for the interactions listed in [pairs]
> >>> that do not have corresponding [pairtypes]. FudgeLJ and
> >>> [nonbond_params] are used in such generation, per other parts of 5.7.
> >>>
>    Or am I completely wrong?
>  In my set up then, are 1-4 Coulombic interactions determined by the
>  pair
>  list and fudge QQ?
> >>> If the contradiction you think exists is this one...
> > On 02/08/11, *Gavin Melaugh *  wrote:
> >> Hi Justin
> >>
> >> Again thanks for the reply. I am not disagreeing with you but If I
>  don't
> >> include a [pairs] directive in the topology file (with gen_pairs
>  =yes),
> >> then there are no 1-4 LJ nor 1-4 Coulombic energies written in
> >> the log
> >> file. When I include the [pair s] directive then both types of
> >> interaction are written to the log file. Therefore does gen_pairs=
>  yes +
> >> [pairs] directive generate 1,4 LJ and 1,4 Coulomb according to
>  fudge LJ
> >> and QQ?
> >>> ... then 5.3.4 indicates that the presence of a [pairs] directive will
> >>> generate the 1,4 output fields. The parameters for that output are
> >>> taken from [pairtypes]. If gen-pairs=yes then the parameters are
> >>> generated, else some warning/error occurs. The example in 5.7.1 has
> >

Re: [gmx-users] atom types

[Gavin Melaugh]

I hope this doesn't come across as stupid, or worse insolent. But what
is the point in stating the charges in the atom type section then ?

Gavin

Mark Abraham wrote:
>
>
> On 02/08/11, *Gavin Melaugh *  wrote:
>> Hi Justin
>>
>> I see that we may have got our wires crossed from the off.
>> Consider the [pairs] directive, which determines which atoms interact
>> in a 1,4 manner.  Consider two atoms listed the [pairs] directive. From
>> the point of the Coulombic interaction between these two atoms I suppose
>> my original question should have been: Does mdrun, when calculating the
>> Coulombic potential between these two atoms, use the charges assigned to
>> the atoms in the [atomtypes] directive or [atoms] directive ?
>
> The charges from [atoms] are used. The charges in [atomtypes] are
> never used in any forcefield currently used with GROMACS. It might not
> be possible for grompp to ever use them, but I'd have to check the
> code for that.
>
> Mark
>
>>
>>
>> Cheers
>>
>> Gavin
>>  
>> Justin A. Lemkul wrote:
>> >
>> >
>> > Gavin Melaugh wrote:
>> >> Yes I think the example vindicates what I am saying as well. I
>> suppose I
>> >> the "contradiction" ( I'll call it the point of confusion) you
>> refer to
>> >> is perhaps when Justin (who is always more than helpful) said that
>> >> "Charges are irrelevant for generation of pair interactions. 
>> Nonbonded
>> >> pair interactions are LJ, not Coulombic.  You will certainly have 1-4
>> >> Coulombic interactions, but they are not generated by gen_pairs.  See
>> >> manual section 5.3.4."
>> >>
>> >
>> > Charges *are* irrelevant - the information is not used when generating
>> > pairs, which I thought was the original question.  The charge
>> > information is used during MD, when the pair list tells mdrun which
>> > atoms interact in what way. Then you get Coul. 1-4 terms.  Perhaps I
>> > missed your point, but this whole thread started as "which charge is
>> > used to generate pairs?"  The answer is still none.
>> >
>> >> My sequence of 1-4 interaction generation should go like this I
>> suppose:
>> >>
>> >> e.g.
>> >> [pairs]
>> >> 3  6
>> >> no parameters present  therefore get from [pairtypes] directive.
>> >> no [pairtypes] directive therefore get from [non_bonded parameters]
>> >> directive as gen pairs = yes
>> >> again no [non_bonded parameters] directive.
>> >> Therefore generate 1,4 interaction parameters based on the normal
>> sigma
>> >> and epsilon values (comb rule 3) present in [atomtypes] directive, in
>> >> accordance with fudge LJ and QQ.
>> >>
>> >> My point is, then in conclusion, that in this way surely the 1,4
>> >> electrostatic interactions are determined by the pair list and in my
>> >> case gen_pairs = yes no?
>> >>
>> >
>> > Yes.
>> >
>> > -Justin
>> >
>> >> Many Thanks
>> >>
>> >> Gavin
>> >>
>> >>
>> >> Mark Abraham wrote:
>> >>> On 02/08/11, *Gavin Melaugh *  wrote:
>>  Hi Mark
>> 
>> 
>>  Thanks for the reply.
>>  I am currently reading that section of the manual and, unless I am
>>  completely mistaken, it seems to vindicate what I am saying.
>>  "Extra Lennard-Jones and electrostatic interactions between pairs of
>>  atoms in a molecule can be added in the [pairs] section of a
>> molecule
>>  definition".
>>  In my [atom types] directive I have atomtype, charge mass, sigma and
>>  epsilon etc. All nonbonding parameters are then calculated
>>  according to
>>  the combination rule (in my case 3). 1-4 interactions are then
>>  calculated based on the information in [pairs] directive (all atoms
>>  are
>>  three bond away). I just have the atom indices of each pair in this
>>  directive therefore with gen_pairs = yes, the interaction parameters
>>  between each pair (which are 1-4) are calculated based on Fudge
>> LJ and
>>  Fudge QQ (which are both 0.5 in my case). All of this in conjunction
>>  with nrexcl =3.
>> >>> That will generate parameters for the interactions listed in [pairs]
>> >>> that do not have corresponding [pairtypes]. FudgeLJ and
>> >>> [nonbond_params] are used in such generation, per other parts of 5.7.
>> >>>
>>    Or am I completely wrong?
>>  In my set up then, are 1-4 Coulombic interactions determined by the
>>  pair
>>  list and fudge QQ?
>> >>> If the contradiction you think exists is this one...
>> > On 02/08/11, *Gavin Melaugh *  wrote:
>> >> Hi Justin
>> >>
>> >> Again thanks for the reply. I am not disagreeing with you but If I
>>  don't
>> >> include a [pairs] directive in the topology file (with gen_pairs
>>  =yes),
>> >> then there are no 1-4 LJ nor 1-4 Coulombic energies written in
>> >> the log
>> >> file. When I include the [pair s] directive then both types of
>> >> interaction are written to the log file. Therefore does gen_pairs=
>>  yes +
>> >> [pairs] directive generate 1,4 LJ and 1,4 Coulomb according to
>>  fudge LJ
>> >> an

[gmx-users] pair types in FEP calculations

[Marcelino Arciniega Castro]

Hi there,

Could someone please help to understand the use of function types 3 and 2 in 
the 
[pairtypes ] directive in the topology file?
I am trying to used it in a free energy calculation but the gromacs 
manual-4.5.4 in the
page 112 is kind of confusing.
Thanks in advance,
Marcelino 
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] atom types

[Justin A. Lemkul]

Gavin Melaugh wrote:

I hope this doesn't come across as stupid, or worse insolent. But what
is the point in stating the charges in the atom type section then ?



They're probably a holdover from some earlier Gromacs version (an ancient one) 
that used those charges somehow, or they were listed there as some grand idea to 
have a more universal, streamlined topology builder that could use generic 
charges.  Just a guess.  They stayed there because it's far easier to just 
ignore the unnecessary information than the re-write the code and those files 
such that the lines are read properly with one less field.  The charges can be 
useful for starting parameterization; they are based on common functional groups 
that may not be far off from the correct values in many cases.


-Justin


Gavin

Mark Abraham wrote:


On 02/08/11, *Gavin Melaugh *  wrote:

Hi Justin

I see that we may have got our wires crossed from the off.
Consider the [pairs] directive, which determines which atoms interact
in a 1,4 manner.  Consider two atoms listed the [pairs] directive. From
the point of the Coulombic interaction between these two atoms I suppose
my original question should have been: Does mdrun, when calculating the
Coulombic potential between these two atoms, use the charges assigned to
the atoms in the [atomtypes] directive or [atoms] directive ?

The charges from [atoms] are used. The charges in [atomtypes] are
never used in any forcefield currently used with GROMACS. It might not
be possible for grompp to ever use them, but I'd have to check the
code for that.

Mark



Cheers

Gavin
 
Justin A. Lemkul wrote:


Gavin Melaugh wrote:

Yes I think the example vindicates what I am saying as well. I

suppose I

the "contradiction" ( I'll call it the point of confusion) you

refer to

is perhaps when Justin (who is always more than helpful) said that
"Charges are irrelevant for generation of pair interactions. 

Nonbonded

pair interactions are LJ, not Coulombic.  You will certainly have 1-4
Coulombic interactions, but they are not generated by gen_pairs.  See
manual section 5.3.4."


Charges *are* irrelevant - the information is not used when generating
pairs, which I thought was the original question.  The charge
information is used during MD, when the pair list tells mdrun which
atoms interact in what way. Then you get Coul. 1-4 terms.  Perhaps I
missed your point, but this whole thread started as "which charge is
used to generate pairs?"  The answer is still none.


My sequence of 1-4 interaction generation should go like this I

suppose:

e.g.
[pairs]
3  6
no parameters present  therefore get from [pairtypes] directive.
no [pairtypes] directive therefore get from [non_bonded parameters]
directive as gen pairs = yes
again no [non_bonded parameters] directive.
Therefore generate 1,4 interaction parameters based on the normal

sigma

and epsilon values (comb rule 3) present in [atomtypes] directive, in
accordance with fudge LJ and QQ.

My point is, then in conclusion, that in this way surely the 1,4
electrostatic interactions are determined by the pair list and in my
case gen_pairs = yes no?


Yes.

-Justin


Many Thanks

Gavin


Mark Abraham wrote:

On 02/08/11, *Gavin Melaugh *  wrote:

Hi Mark


Thanks for the reply.
I am currently reading that section of the manual and, unless I am
completely mistaken, it seems to vindicate what I am saying.
"Extra Lennard-Jones and electrostatic interactions between pairs of
atoms in a molecule can be added in the [pairs] section of a

molecule

definition".
In my [atom types] directive I have atomtype, charge mass, sigma and
epsilon etc. All nonbonding parameters are then calculated
according to
the combination rule (in my case 3). 1-4 interactions are then
calculated based on the information in [pairs] directive (all atoms
are
three bond away). I just have the atom indices of each pair in this
directive therefore with gen_pairs = yes, the interaction parameters
between each pair (which are 1-4) are calculated based on Fudge

LJ and

Fudge QQ (which are both 0.5 in my case). All of this in conjunction
with nrexcl =3.

That will generate parameters for the interactions listed in [pairs]
that do not have corresponding [pairtypes]. FudgeLJ and
[nonbond_params] are used in such generation, per other parts of 5.7.


  Or am I completely wrong?
In my set up then, are 1-4 Coulombic interactions determined by the
pair
list and fudge QQ?

If the contradiction you think exists is this one...

On 02/08/11, *Gavin Melaugh *  wrote:

Hi Justin

Again thanks for the reply. I am not disagreeing with you but If I

don't

include a [pairs] directive in the topology file (with gen_pairs

=yes),

then there are no 1-4 LJ nor 1-4 Coulombic energies written in
the log
file. When I include the [pair s] directive then both types of
interaction are written to the log file. Therefore does gen_pairs=

yes +

[pairs] directive generate 1,4 LJ and 1,4 Coulomb according to

fudge LJ

and QQ?

... then 

Re: [gmx-users] pair types in FEP calculations

[Justin A. Lemkul]

Marcelino Arciniega Castro wrote:

Hi there,

Could someone please help to understand the use of function types 3 and 2 in the 
[pairtypes ] directive in the topology file?

I am trying to used it in a free energy calculation but the gromacs 
manual-4.5.4 in the
page 112 is kind of confusing.


That information may be outdated.  The .mdp keyword "couple-intramol" seems to 
have the same function, negating the need to change anything in the [pairtypes]. 
 Just a guess.  Please test and see if your results are reasonable.


-Justin

--


Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.

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

Re: [gmx-users] atom types

[Mark Abraham]

On 02/08/11, Gavin Melaugh   wrote:
> I hope this doesn't come across as stupid, or worse insolent. But what
> is the point in stating the charges in the atom type section then ?
> 

It would make sense for a force field that had a fixed charge for some/all atom 
types. If any such force fields exist, they don't get used for biomolecular MD. 
I don't know if such a force field can be implemented in GROMACS (per my 
earlier comment about checking the code). In some cases, grompp does a bunch of 
fancy footwork to try to infer the format of the content of the input lines. 
It's possible it can recognize that a charge field is missing and perhaps the 
charge gets read from [atomtypes] - but that is just not useful for 
biomolecular MD.

Another possibility is that the file formats serve other purposes (possibly in 
the distant past).

Mark


> 
> 
> Gavin
> 
> Mark Abraham wrote:
> >
> >
> > On 02/08/11, *Gavin Melaugh *  wrote:
> >> Hi Justin
> >>
> >> I see that we may have got our wires crossed from the off.
> >> Consider the [pairs] directive, which determines which atoms interact
> >> in a 1,4 manner.  Consider two atoms listed the [pairs] directive. From
> >> the point of the Coulombic interaction between these two atoms I suppose
> >> my original question should have been: Does mdrun, when calculating the
> >> Coulombic potential between these two atoms, use the charges assigned to
> >> the atoms in the [atomtypes] directive or [atoms] directive ?
> >
> > The charges from [atoms] are used. The charges in [atomtypes] are
> > never used in any forcefield currently used with GROMACS. It might not
> > be possible for grompp to ever use them, but I'd have to check the
> > code for that.
> >
> > Mark
> >
> >>
> >>
> >> Cheers
> >>
> >> Gavin
> >>  
> >> Justin A. Lemkul wrote:
> >> >
> >> >
> >> > Gavin Melaugh wrote:
> >> >> Yes I think the example vindicates what I am saying as well. I
> >> suppose I
> >> >> the "contradiction" ( I'll call it the point of confusion) you
> >> refer to
> >> >> is perhaps when Justin (who is always more than helpful) said that
> >> >> "Charges are irrelevant for generation of pair interactions. 
> >> Nonbonded
> >> >> pair interactions are LJ, not Coulombic.  You will certainly have 1-4
> >> >> Coulombic interactions, but they are not generated by gen_pairs.  See
> >> >> manual section 5.3.4."
> >> >>
> >> >
> >> > Charges *are* irrelevant - the information is not used when generating
> >> > pairs, which I thought was the original question.  The charge
> >> > information is used during MD, when the pair list tells mdrun which
> >> > atoms interact in what way. Then you get Coul. 1-4 terms.  Perhaps I
> >> > missed your point, but this whole thread started as "which charge is
> >> > used to generate pairs?"  The answer is still none.
> >> >
> >> >> My sequence of 1-4 interaction generation should go like this I
> >> suppose:
> >> >>
> >> >> e.g.
> >> >> [pairs]
> >> >> 3  6
> >> >> no parameters present  therefore get from [pairtypes] directive.
> >> >> no [pairtypes] directive therefore get from [non_bonded parameters]
> >> >> directive as gen pairs = yes
> >> >> again no [non_bonded parameters] directive.
> >> >> Therefore generate 1,4 interaction parameters based on the normal
> >> sigma
> >> >> and epsilon values (comb rule 3) present in [atomtypes] directive, in
> >> >> accordance with fudge LJ and QQ.
> >> >>
> >> >> My point is, then in conclusion, that in this way surely the 1,4
> >> >> electrostatic interactions are determined by the pair list and in my
> >> >> case gen_pairs = yes no?
> >> >>
> >> >
> >> > Yes.
> >> >
> >> > -Justin
> >> >
> >> >> Many Thanks
> >> >>
> >> >> Gavin
> >> >>
> >> >>
> >> >> Mark Abraham wrote:
> >> >>> On 02/08/11, *Gavin Melaugh *  wrote:
> >>  Hi Mark
> >> 
> >> 
> >>  Thanks for the reply.
> >>  I am currently reading that section of the manual and, unless I am
> >>  completely mistaken, it seems to vindicate what I am saying.
> >>  "Extra Lennard-Jones and electrostatic interactions between pairs of
> >>  atoms in a molecule can be added in the [pairs] section of a
> >> molecule
> >>  definition".
> >>  In my [atom types] directive I have atomtype, charge mass, sigma and
> >>  epsilon etc. All nonbonding parameters are then calculated
> >>  according to
> >>  the combination rule (in my case 3). 1-4 interactions are then
> >>  calculated based on the information in [pairs] directive (all atoms
> >>  are
> >>  three bond away). I just have the atom indices of each pair in this
> >>  directive therefore with gen_pairs = yes, the interaction parameters
> >>  between each pair (which are 1-4) are calculated based on Fudge
> >> LJ and
> >>  Fudge QQ (which are both 0.5 in my case). All of this in conjunction
> >>  with nrexcl =3.
> >> >>> That will generate parameters for the interactions listed in [pairs]
> >> >>> that do not have corresponding [pa

Re: [gmx-users] pair types in FEP calculations

[Mark Abraham]

On 02/08/11, Marcelino Arciniega Castro   wrote:
> Hi there,
> 
> Could someone please help to understand the use of function types 3 and 2 in 
> the 
> [pairtypes ] directive in the topology file?
> I am trying to used it in a free energy calculation but the gromacs 
> manual-4.5.4 in the
> page 112 is kind of confusing.
> 

That text mentions a bunch of sub-cases. Without knowing what you want to 
achieve and what you find unclear, it's hard to offer any help.

Mark
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

[gmx-users] center of box protein

[nahren manuel]

Dear GMX users,

I am trying to center my protein at the center of box using editconf. When I 
view the same using ngmx, the protein seems to lie outside the box

neweditconf -f alapdb.pdb -bt cubic  -o boxpdb.pdb -c -center -0.156 -0.061 
0.084 -d 1.5

this one works, but the molecule is shifted.

neweditconf -f alapdb.pdb -bt cubic -d 1.5 -o try.pdb

I want my molecule (geometrical )  center to be the center of the box with -d 
1.5.


Best,
nahren
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] center of box protein

[Justin A. Lemkul]

nahren manuel wrote:

Dear GMX users,

I am trying to center my protein at the center of box using editconf. 
When I view the same using ngmx, the protein seems to lie outside the box


neweditconf -f alapdb.pdb -bt cubic  -o boxpdb.pdb -c -center -0.156 
-0.061 0.084 -d 1.5




The protein appears "outside" of the box (which is irrelevant in a periodic 
system, anyway) because you're putting it there.  You've specified negative 
coordinates.  Gromacs builds all boxes from the coordinate origin into the 
quadrant where (x,y,z) are all positive.



this one works, but the molecule is shifted.
neweditconf -f alapdb.pdb -bt cubic -d 1.5 -o try.pdb

I want my molecule (geometrical )  center to be the center of the box 
with -d 1.5.





Just use editconf -c -d 1.5; I don't know what "shifted" means above, because 
the combination of -c and -d should center the protein in the box (in fact, the 
use of -d automatically triggers -c).  Using -center overrides it with whatever 
you specify and is only necessary for custom placement.


-Justin


Best,
nahren



--


Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.

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

Re: [gmx-users] malloc on mac

[Sara baretller]

I have a question about the g-hbond selection . after i typed in the g_hbond
, it asked me to choose two groops and i am wondering what is difference
between Protein and Protein-H ?

Thank you

Specify 2 groups to analyze:
Group 0 ( System) has  7722 elements
Group 1 (Protein) has  2000 elements
Group 2 (  Protein-H) has  2000 elements
Group 3 (C-alpha) has 0 elements
Group 4 (   Backbone) has 0 elements
Group 5 (  MainChain) has 0 elements
Group 6 (   MainChain+Cb) has 0 elements
Group 7 (MainChain+H) has 0 elements
Group 8 (  SideChain) has  2000 elements
Group 9 (SideChain-H) has  2000 elements
Group10 (Prot-Masses) has  2000 elements
Group11 (non-Protein) has  5722 elements
Group12 (  Other) has  5722 elements
Group13 (  W) has  5722 elemen




On Thu, Jul 14, 2011 at 11:10 PM, Mark Abraham wrote:

> On 15/07/2011 11:48 AM, Itamar Kass wrote:
>
>> Hi all,
>>
>> I am trying to find all possible h-bonds between chains in my complex. I
>> am using:
>>
>> g_hbond_d -f system_run1_MD050_fitbb.xtc -s system_for_EM.tpr -n
>> system.ndx -g system_run1_MD050_fitbb_Hbond_**all.log -num
>> system_run1_MD050_fitbb_Hbnum_**all.xvg -hbn system_run1_MD050_fitbb.xtc_
>> **Hbond_all.ndx
>>
>> and after reading 100ns the system crash and report:
>>
>> Reading frame1000 time 10.000   g_hbond_d(42677) malloc: *** error
>> for object 0xd: pointer being reallocated was not allocated
>> *** set a breakpoint in malloc_error_break to debug
>>
>> I am suing gromacs 4.0.7 (double precision) on mac (10.6.8), which I
>> compiled myself (not via fink). I also tried to compile it using dmalloc
>> without success.
>>
>
> Sounds like a bug. Try g_hbond from a more recent version of GROMACS.
>
> BTW, you should only compile and use tools in double precision if that
> extra accuracy is what you actually need. Otherwise you're just running
> slower than you need to. All tools can read files written at either
> precision.
>
> Mark
>
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/**mailman/listinfo/gmx-users
> Please search the archive at http://www.gromacs.org/**
> Support/Mailing_Lists/Searchbefore
>  posting!
> Please don't post (un)subscribe requests to the list. Use the www interface
> or send it to gmx-users-requ...@gromacs.org.
> Can't post? Read 
> http://www.gromacs.org/**Support/Mailing_Lists
>
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] malloc on mac

[Justin A. Lemkul]

Sara baretller wrote:
I have a question about the g-hbond selection . after i typed in the 
g_hbond , it asked me to choose two groops and i am wondering what is 
difference between Protein and Protein-H ?




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

-Justin


Thank you

Specify 2 groups to analyze:
Group 0 ( System) has  7722 elements
Group 1 (Protein) has  2000 elements
Group 2 (  Protein-H) has  2000 elements
Group 3 (C-alpha) has 0 elements
Group 4 (   Backbone) has 0 elements
Group 5 (  MainChain) has 0 elements
Group 6 (   MainChain+Cb) has 0 elements
Group 7 (MainChain+H) has 0 elements
Group 8 (  SideChain) has  2000 elements
Group 9 (SideChain-H) has  2000 elements
Group10 (Prot-Masses) has  2000 elements
Group11 (non-Protein) has  5722 elements
Group12 (  Other) has  5722 elements
Group13 (  W) has  5722 elemen




On Thu, Jul 14, 2011 at 11:10 PM, Mark Abraham > wrote:


On 15/07/2011 11:48 AM, Itamar Kass wrote:

Hi all,

I am trying to find all possible h-bonds between chains in my
complex. I am using:

g_hbond_d -f system_run1_MD050_fitbb.xtc -s system_for_EM.tpr -n
system.ndx -g system_run1_MD050_fitbb_Hbond___all.log -num
system_run1_MD050_fitbb_Hbnum___all.xvg -hbn
system_run1_MD050_fitbb.xtc___Hbond_all.ndx

and after reading 100ns the system crash and report:

Reading frame1000 time 10.000   g_hbond_d(42677) malloc:
*** error for object 0xd: pointer being reallocated was not
allocated
*** set a breakpoint in malloc_error_break to debug

I am suing gromacs 4.0.7 (double precision) on mac (10.6.8),
which I compiled myself (not via fink). I also tried to compile
it using dmalloc without success.


Sounds like a bug. Try g_hbond from a more recent version of GROMACS.

BTW, you should only compile and use tools in double precision if
that extra accuracy is what you actually need. Otherwise you're just
running slower than you need to. All tools can read files written at
either precision.

Mark

-- 
gmx-users mailing listgmx-users@gromacs.org


http://lists.gromacs.org/__mailman/listinfo/gmx-users

Please search the archive at
http://www.gromacs.org/__Support/Mailing_Lists/Search
 before posting!
Please don't post (un)subscribe requests to the list. Use the www
interface or send it to gmx-users-requ...@gromacs.org
.
Can't post? Read http://www.gromacs.org/__Support/Mailing_Lists





--


Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.

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

[gmx-users] Heat of Vaporization

[Elisabeth]

Hello,

I wanted to know your ideas on calculation of heat of vaporization using a
single phase run rather than running two separate simulations for liquid and
gas!

1- Two separate simulations for liquid and gas

DHvap =  -  + RT

1a:  - <*total* potential of a single chain in vacu
> ( bond+angle+torsion + nonbonded interaction of chain with itself)

or

 1b: * - < intra potential of a single chain in vacu
> ( bond+angle+torsion) *

2- Single liquid phase run: (non need to run in vacu)

2a : DHvap =  - < intra molecular potential terms in
liquid phase> (same liquid phase simulation by adding up bond+angle+torsion
terms)

2a: In other worlds *DHvap=  *

In my case the latter definition is giving much more accurate results than
1a.

I would like to know your idea and comments on methods 1b and 2a.

Appreciate your comments.

Best,
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] malloc on mac

[Sara baretller]

does the g_hbond work for course grained file. i tried this command and it
gave me nothing but an error it does not matter wich selection



On Tue, Aug 2, 2011 at 1:01 PM, Justin A. Lemkul  wrote:

>
>
> Sara baretller wrote:
>
>> I have a question about the g-hbond selection . after i typed in the
>> g_hbond , it asked me to choose two groops and i am wondering what is
>> difference between Protein and Protein-H ?
>>
>>
> http://www.gromacs.org/**Documentation/Terminology/**Default_Index_Groups
>
> -Justin
>
>  Thank you
>>
>> Specify 2 groups to analyze:
>> Group 0 ( System) has  7722 elements
>> Group 1 (Protein) has  2000 elements
>> Group 2 (  Protein-H) has  2000 elements
>> Group 3 (C-alpha) has 0 elements
>> Group 4 (   Backbone) has 0 elements
>> Group 5 (  MainChain) has 0 elements
>> Group 6 (   MainChain+Cb) has 0 elements
>> Group 7 (MainChain+H) has 0 elements
>> Group 8 (  SideChain) has  2000 elements
>> Group 9 (SideChain-H) has  2000 elements
>> Group10 (Prot-Masses) has  2000 elements
>> Group11 (non-Protein) has  5722 elements
>> Group12 (  Other) has  5722 elements
>> Group13 (  W) has  5722 elemen
>>
>>
>>
>>
>> On Thu, Jul 14, 2011 at 11:10 PM, Mark Abraham 
>> > mark.abra...@anu.edu.**au >> wrote:
>>
>>On 15/07/2011 11:48 AM, Itamar Kass wrote:
>>
>>Hi all,
>>
>>I am trying to find all possible h-bonds between chains in my
>>complex. I am using:
>>
>>g_hbond_d -f system_run1_MD050_fitbb.xtc -s system_for_EM.tpr -n
>>system.ndx -g system_run1_MD050_fitbb_Hbond_**__all.log -num
>>system_run1_MD050_fitbb_Hbnum_**__all.xvg -hbn
>>system_run1_MD050_fitbb.xtc___**Hbond_all.ndx
>>
>>and after reading 100ns the system crash and report:
>>
>>Reading frame1000 time 10.000   g_hbond_d(42677) malloc:
>>*** error for object 0xd: pointer being reallocated was not
>>allocated
>>*** set a breakpoint in malloc_error_break to debug
>>
>>I am suing gromacs 4.0.7 (double precision) on mac (10.6.8),
>>which I compiled myself (not via fink). I also tried to compile
>>it using dmalloc without success.
>>
>>
>>Sounds like a bug. Try g_hbond from a more recent version of GROMACS.
>>
>>BTW, you should only compile and use tools in double precision if
>>that extra accuracy is what you actually need. Otherwise you're just
>>running slower than you need to. All tools can read files written at
>>either precision.
>>
>>Mark
>>
>>-- gmx-users mailing listgmx-users@gromacs.org
>>
>>
>>
>> http://lists.gromacs.org/__**mailman/listinfo/gmx-users
>>
>> 
>> >
>>Please search the archive at
>>
>> http://www.gromacs.org/__**Support/Mailing_Lists/Search
>>
>> >
>> before posting!
>>Please don't post (un)subscribe requests to the list. Use the www
>>interface or send it to gmx-users-requ...@gromacs.org
>>> >.
>>
>>Can't post? Read 
>> http://www.gromacs.org/__**Support/Mailing_Lists
>>
>> 
>> >
>>
>>
>>
> --
> ==**==
>
> Justin A. Lemkul
> Ph.D. Candidate
> ICTAS Doctoral Scholar
> MILES-IGERT Trainee
> Department of Biochemistry
> Virginia Tech
> Blacksburg, VA
> jalemkul[at]vt.edu | (540) 231-9080
> http://www.bevanlab.biochem.**vt.edu/Pages/Personal/justin
>
> ==**==
>
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/**mailman/listinfo/gmx-users
> Please search the archive at http://www.gromacs.org/**
> Support/Mailing_Lists/Searchbefore
>  posting!
> Please don't post (un)subscribe requests to the list. Use the www interface
> or send it to gmx-users-requ...@gromacs.org.
> Can't post? Read 
> http://www.gromacs.org/**Support/Mailing_Lists
>
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post 

Re: [gmx-users] malloc on mac

[Justin A. Lemkul]

Sara baretller wrote:


does the g_hbond work for course grained file. i tried this command and 
it gave me nothing but an error it does not matter wich selection





If there are no hydrogens, there is nothing for g_hbond to measure.  You need an 
atomistic representation such that D-H-A angles and distances can be measured.


-Justin



On Tue, Aug 2, 2011 at 1:01 PM, Justin A. Lemkul > wrote:




Sara baretller wrote:

I have a question about the g-hbond selection . after i typed in
the g_hbond , it asked me to choose two groops and i am
wondering what is difference between Protein and Protein-H ?


http://www.gromacs.org/__Documentation/Terminology/__Default_Index_Groups


-Justin

Thank you

Specify 2 groups to analyze:
Group 0 ( System) has  7722 elements
Group 1 (Protein) has  2000 elements
Group 2 (  Protein-H) has  2000 elements
Group 3 (C-alpha) has 0 elements
Group 4 (   Backbone) has 0 elements
Group 5 (  MainChain) has 0 elements
Group 6 (   MainChain+Cb) has 0 elements
Group 7 (MainChain+H) has 0 elements
Group 8 (  SideChain) has  2000 elements
Group 9 (SideChain-H) has  2000 elements
Group10 (Prot-Masses) has  2000 elements
Group11 (non-Protein) has  5722 elements
Group12 (  Other) has  5722 elements
Group13 (  W) has  5722 elemen




On Thu, Jul 14, 2011 at 11:10 PM, Mark Abraham
mailto:mark.abra...@anu.edu.au>
>> wrote:

   On 15/07/2011 11:48 AM, Itamar Kass wrote:

   Hi all,

   I am trying to find all possible h-bonds between chains in my
   complex. I am using:

   g_hbond_d -f system_run1_MD050_fitbb.xtc -s
system_for_EM.tpr -n
   system.ndx -g system_run1_MD050_fitbb_Hbond_all.log -num
   system_run1_MD050_fitbb_Hbnum_all.xvg -hbn
   system_run1_MD050_fitbb.xtc_Hbond_all.ndx

   and after reading 100ns the system crash and report:

   Reading frame1000 time 10.000   g_hbond_d(42677)
malloc:
   *** error for object 0xd: pointer being reallocated was not
   allocated
   *** set a breakpoint in malloc_error_break to debug

   I am suing gromacs 4.0.7 (double precision) on mac (10.6.8),
   which I compiled myself (not via fink). I also tried to
compile
   it using dmalloc without success.


   Sounds like a bug. Try g_hbond from a more recent version of
GROMACS.

   BTW, you should only compile and use tools in double precision if
   that extra accuracy is what you actually need. Otherwise
you're just
   running slower than you need to. All tools can read files
written at
   either precision.

   Mark

   -- gmx-users mailing listgmx-users@gromacs.org

   >

   http://lists.gromacs.org/mailman/listinfo/gmx-users

   >
   Please search the archive at
   http://www.gromacs.org/Support/Mailing_Lists/Search

   > before
posting!
   Please don't post (un)subscribe requests to the list. Use the www
   interface or send it to gmx-users-requ...@gromacs.org

   >.

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

   >



-- 
==__==


Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu  | (540) 231-9080

http://www.bevanlab.biochem.__vt.edu/Pages/Personal/justin

Re: [gmx-users] malloc on mac

[Tsjerk Wassenaar]

Hi Sara,

Since you don't have hydrogens, it's hard to talk about hydrogen bonds
according to conventional criteria.

Note this has nothing to do with malloc; please start a new thread rather
than asking your question in another.

Cheers,

Tsjerk

On Aug 2, 2011 7:56 PM, "Sara baretller" 
wrote:


does the g_hbond work for course grained file. i tried this command and it
gave me nothing but an error it does not matter wich selection

On Tue, Aug 2, 2011 at 1:01 PM, Justin A. Lemkul  wrote: >
> > > Sara baretller...

--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] malloc on mac

[Sara baretller]

do you suggest any other way to find the hbond?? is the cut of distance of 3
or less something i can do in this case or not??

thank you

On Tue, Aug 2, 2011 at 1:58 PM, Justin A. Lemkul  wrote:

>
>
> Sara baretller wrote:
>
>>
>> does the g_hbond work for course grained file. i tried this command and it
>> gave me nothing but an error it does not matter wich selection
>>
>>
>>
> If there are no hydrogens, there is nothing for g_hbond to measure.  You
> need an atomistic representation such that D-H-A angles and distances can be
> measured.
>
> -Justin
>
>
>> On Tue, Aug 2, 2011 at 1:01 PM, Justin A. Lemkul > jalem...@vt.edu>> wrote:
>>
>>
>>
>>Sara baretller wrote:
>>
>>I have a question about the g-hbond selection . after i typed in
>>the g_hbond , it asked me to choose two groops and i am
>>wondering what is difference between Protein and Protein-H ?
>>
>>
>>http://www.gromacs.org/__**Documentation/Terminology/__**
>> Default_Index_Groups
>>> Default_Index_Groups
>> >
>>
>>-Justin
>>
>>Thank you
>>
>>Specify 2 groups to analyze:
>>Group 0 ( System) has  7722 elements
>>Group 1 (Protein) has  2000 elements
>>Group 2 (  Protein-H) has  2000 elements
>>Group 3 (C-alpha) has 0 elements
>>Group 4 (   Backbone) has 0 elements
>>Group 5 (  MainChain) has 0 elements
>>Group 6 (   MainChain+Cb) has 0 elements
>>Group 7 (MainChain+H) has 0 elements
>>Group 8 (  SideChain) has  2000 elements
>>Group 9 (SideChain-H) has  2000 elements
>>Group10 (Prot-Masses) has  2000 elements
>>Group11 (non-Protein) has  5722 elements
>>Group12 (  Other) has  5722 elements
>>Group13 (  W) has  5722 elemen
>>
>>
>>
>>
>>On Thu, Jul 14, 2011 at 11:10 PM, Mark Abraham
>>> > >
>>>>>
>> wrote:
>>
>>   On 15/07/2011 11:48 AM, Itamar Kass wrote:
>>
>>   Hi all,
>>
>>   I am trying to find all possible h-bonds between chains in
>> my
>>   complex. I am using:
>>
>>   g_hbond_d -f system_run1_MD050_fitbb.xtc -s
>>system_for_EM.tpr -n
>>   system.ndx -g system_run1_MD050_fitbb_Hbond_**all.log
>> -num
>>   system_run1_MD050_fitbb_Hbnum_**all.xvg -hbn
>>   system_run1_MD050_fitbb.xtc___**__Hbond_all.ndx
>>
>>   and after reading 100ns the system crash and report:
>>
>>   Reading frame1000 time 10.000   g_hbond_d(42677)
>>malloc:
>>   *** error for object 0xd: pointer being reallocated was not
>>   allocated
>>   *** set a breakpoint in malloc_error_break to debug
>>
>>   I am suing gromacs 4.0.7 (double precision) on mac (10.6.8),
>>   which I compiled myself (not via fink). I also tried to
>>compile
>>   it using dmalloc without success.
>>
>>
>>   Sounds like a bug. Try g_hbond from a more recent version of
>>GROMACS.
>>
>>   BTW, you should only compile and use tools in double precision
>> if
>>   that extra accuracy is what you actually need. Otherwise
>>you're just
>>   running slower than you need to. All tools can read files
>>written at
>>   either precision.
>>
>>   Mark
>>
>>   -- gmx-users mailing listgmx-users@gromacs.org
>>
>>   **>
>>
>>
>>   
>> http://lists.gromacs.org/**mailman/listinfo/gmx-users
>>
>> 
>> >
>>   
>> 
>>
>> 
>> >>
>>   Please search the archive at
>>   
>> http://www.gromacs.org/**Support/Mailing_Lists/Search
>>
>> 
>> >
>>   
>> 

[gmx-users] (no subject)

[Sara baretller]

Hi Justin you mentioned in the past email that i can not use the measure
hbond;
do you suggest any other way to find the hbond??

thank you
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

[gmx-users] Re: coarse grain hydrogen bonding

[Justin A. Lemkul]

Tsjerk's message covered this:

http://lists.gromacs.org/pipermail/gmx-users/2011-August/063472.html

To detect a hydrogen bond in a conventional sense using geometric criteria, you 
need an atomistic representation.  There are tools available to reconstruct 
atomic positions from CG models, but I have never used them so I will not 
comment further other than to suggest you may want to investigate them.  Of 
course, that leaves open for debate whether or not the hydrogen bond would 
actually be there or if it is caused solely by the reconstruction algorithm.


My sense is that no, there isn't a very easy way to collect hydrogen bonding 
data from a CG trajectory.


-Justin

Sara baretller wrote:
Hi Justin you mentioned in the past email that i can not use the measure 
hbond;

do you suggest any other way to find the hbond??

thank you


--


Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.

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

[gmx-users] Re: Heat of Vaporization

[Dr. Vitaly V. Chaban]

>
> Hello,
>
> I wanted to know your ideas on calculation of heat of vaporization using a
> single phase run rather than running two separate simulations for liquid and
> gas!
>
> 1- Two separate simulations for liquid and gas
>
> DHvap =  -  + RT
>
> 1a:  - <*total* potential of a single chain in vacu
>> ( bond+angle+torsion + nonbonded interaction of chain with itself)
>
> or
>
>  1b: * - < intra potential of a single chain in vacu
>> ( bond+angle+torsion) *
>
> 2- Single liquid phase run: (non need to run in vacu)
>
> 2a : DHvap =  - < intra molecular potential terms in
> liquid phase> (same liquid phase simulation by adding up bond+angle+torsion
> terms)
>
> 2a: In other worlds *DHvap=  *
>
> In my case the latter definition is giving much more accurate results than
> 1a.
>
> I would like to know your idea and comments on methods 1b and 2a.
>
> Appreciate your comments.
>

What is your system?


(1a) is what everyone does.

-- 
Dr. Vitaly V. Chaban, 430 Hutchison Hall, Chem. Dept.
Univ. Rochester, Rochester, New York 14627-0216
THE UNITED STATES OF AMERICA
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

[gmx-users] Jarzynski and PMF

[Sai Kumar Ramadugu]

Dear All,

I have tried to calculate the free energy of binding using Jarzynski
equality. I employed the following procedure.

--I did force pulling simulations along Z-direction as exemplified in
Justin's umbrella sampling simulations. I did 50, 250 and 500 pulling
simulations to test the convergence and also two different pulling rates
(2nm/ns and 4nm/ns), two different force
  constants (2000 kJ mol^-1 nm^-2 and 4000 kJ mol^-1 nm^-2).

--From the pulling simulation, I get the pullf.xvg and pullx.xvg.

--I use the force from pullf.xvg and calculate the work at each time step by
multiplying with v*dt and calculate average work for each simulation.

--Once I have the work for each simulation, I calculate the exp(-beta*W) for
each simulation and calculate the average of exp(-beta*W).

--Once I have the average of exp(-beta*W), I calculate the free energy delF
= -kB*T ln(exp(-beta*W)).

The question I have is whether my approach for calculating the work is
correct/feasible or not.
Does the pullf.xvg written by mdrun during a pulling simulation contain the
pull force only?
Do I need to add any corrections?
I'm doing NPT simulations, so can I add the pdV correction as suggested in
one of the paper (Macromolecules, 2008, 41 (6), pp 2283–2289) by Berk Hess?
Or did I understand the paper incorrectly?

For a publication level of work, I'll use NVT ensemble. But for now my
biggest concern is calculation of work.

Any suggestions will be greatly helpful.

Thanks for you time,

Sai Ramadugu
Dept of Chemistry,
University of Iowa.
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] bug in demux.pl

[wibke . sudholt]

Dear Email Sender,

Thank you very much for contacting me! Unfortunately, I am not available in the 
office at the moment and cannot respond to your email. I will be able to handle 
your request starting again Thursday, August 4, 2011.

For all questions about CloudBroker and the CloudBroker Platform, please 
contact the company under i...@cloudbroker.com or Nicola Fantini under 
nicola.fant...@cloudbroker.com. If you need to talk to me in person for urgent 
or important issues, please call my mobile phone number.

Kind regards,

Wibke Sudholt


-- 
Dr. Wibke Sudholt
CEO
CloudBroker GmbH
Technoparkstrasse 1
CH-8005 Zurich
Switzerland
Phone: +41 44 633 79 34
Email:  i...@cloudbroker.com
Web:http://www.cloudbroker.com



-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

[gmx-users] bug in demux.pl

[Swarnendu Tripathi]

Hi,

I am not sure if this problem has been solved or discussed before.

I found that when I use the demux.pl script it does not read the actual time
from the .log file. This is for gromacs-4.0.7.

These are few lines from one of the .log file for my replica exchnage
simulation. For example in the second line below the actual time is
"112.0" in ps. However, the demux.pl script reads the time from the
"Replica exchange" phrase as "1.1e+06". Same thing repeats in the
following steps which is not right.
==
Step   Time Lambda
  28000   112.125000.0

   Energies (kJ/mol)
G96Bond   G96AngleProper Dih.  Improper Dih.  LJ-14
1.58180e+031.49412e+037.67393e+023.69239e+02   -2.78302e+02
 Coulomb-14LJ (SR)   Coulomb (SR)  PotentialKinetic En.
   -2.39467e+008.01055e+01   -1.05400e+014.00143e+034.06353e+03
   Total EnergyTemperature Pressure (bar)
8.06496e+034.51898e+020.0e+00

Replica exchange at step 28000 time 1.1e+06
Repl 7 <-> 8  dE =  1.607e-01
Repl ex  01 x  23 x  45 x  67 x  89 x 10   11 x 12   13
x 14   15 x 16   17 x 18   19 x 20   21 x 22   23
Repl pr.79   1.0   1.0   .85   1.0   .98
.98   .90   1.0   1.0   .90

   Step   Time Lambda
  29000   116.125000.0

   Energies (kJ/mol)
G96Bond   G96AngleProper Dih.  Improper Dih.  LJ-14
1.42588e+031.43257e+037.64199e+023.39804e+02   -7.50646e+02
 Coulomb-14LJ (SR)   Coulomb (SR)  PotentialKinetic En.
   -2.96172e+009.62908e+01   -1.95152e+013.28563e+034.16581e+03
   Total EnergyTemperature Pressure (bar)
7.45144e+034.63272e+020.0e+00

Replica exchange at step 29000 time 1.2e+06
Repl 8 <-> 9  dE =  2.023e-01
Repl ex  0 x  12 x  34 x  56 x  789   10 x 11   12 x
13   14 x 15   16 x 17   18 x 19   20 x 21   22 x 23
Repl pr   1.0   .98   .66   1.0   .82   1.0
.92   .93   .79   1.0   .96   1.0



As a result my "replica_index.xvg" looks like this which is not right.

1.19996e+06   119   22   12   10   215   23   177   198
31   150   136   20   1442   16   18
1.19996e+069   11   12   22   21   10   2357   178   19
130   156   13   14   2024   18   16
1.19997e+069   12   11   21   22   23   10758   171
19036   15   14   132   20   184   16
1.19997e+06   129   21   11   23   227   10851   17
0   1963   14   152   13   20   18   164
1.19998e+06   12   219   23   11   2278   10510
176   19   1432   15   13   20   16   184
1.19998e+06   21   12   239   22   1187   10501
6   17   19   1423   13   15   16   204   18
1.19998e+06   21   23   12   2298   11   107506
1   19   172   14   133   16   154   20   18
1.1e+06   23   21   22   1298   10   115760
1912   17   13   14   1634   15   18   20
1.1e+06   23   22   219   12   1085   1167   19
021   13   17   16   1443   18   15   20
1.2e+0622   239   21   10   1258   116   1972
0   131   16   174   14   183   20   15
=

Due to the problem in precision, time step repeats. The demux.pl script
should read the actual time from the previous step in the .log file.

-Swarnendu
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

[gmx-users] spatial distribution function in a binary solvent mixture

[dbiswal]

  Dear all,

 I'm working with a binary solvent mixture containing 2000 molecules (1800
 type-A + 200 type-B). Both the types of solvent molecules have similar
structure (they are both diatomic molecules) except the polarity.I'm
trying to calculate sdf of type-B solvent molecules. I followed the
step-by-step instructions from the manual using g_spatial.
 I'm trying to reduce the bin width to 0.05A. The *.cube file generated
with a bin width of 0.09A is already 4.9GB in size. As I'm more
interested in the first solvation shell around the type-B solvent
molecule, I was wondering if I could find a way to control the maximum
radius of the sphere around central molecule (center of coordinates?)?
Also, can anyone please let me know how g_spatial deals with the angular
part of sdf?

(I've searched a lot but could not gather enough meaningful information)

regards,
Debasmita


-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

[gmx-users] spatial distribution function in a binary solvent mixture

[chris . neale]

I think that you have a misconception about what g_spatial does. For a  
system with many type A and many type B, you need to average over all  
of one type as the central solute to compute an rdf, and perhaps that  
is what you want for your sdf. g_spatial, however, does not do any  
fitting. g_sdf did that. You can obtain an old version of the source  
(4.0.7 should have it I think) if you want to use g_sdf. I have never  
used g_sdf myself.


In case that doesn't answer your question, then let me make one more point:

g_spatial requires that a bin exists for every count. Thus either (a)  
find a large memory node somewhere or (b) pre-process your trajectory  
using trjorder to order the solvent molecules and then keep only the N  
closest, then run g_spatial. But again, I think that this will not  
provide what you are looking for but is likely to give you a smear  
over your box.


Chris.

-- original message --

  Dear all,

 I'm working with a binary solvent mixture containing 2000 molecules (1800
 type-A + 200 type-B). Both the types of solvent molecules have similar
structure (they are both diatomic molecules) except the polarity.I'm
trying to calculate sdf of type-B solvent molecules. I followed the
step-by-step instructions from the manual using g_spatial.
 I'm trying to reduce the bin width to 0.05A. The *.cube file generated
with a bin width of 0.09A is already 4.9GB in size. As I'm more
interested in the first solvation shell around the type-B solvent
molecule, I was wondering if I could find a way to control the maximum
radius of the sphere around central molecule (center of coordinates?)?
Also, can anyone please let me know how g_spatial deals with the angular
part of sdf?

(I've searched a lot but could not gather enough meaningful information)

regards,
Debasmita




--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] cannot open file

[Justin A. Lemkul]

Taylor Kaplan wrote:

Hi Justin,

I'm running em.sh file as nohup em.sh &. Below is the code in my 
em.sh file.


mpirun -nolocal -np 8 -machinefile machine mdrun_mpi -nice 0 -v -s 
etExcited.tpr -o etExcited.trr -c etExcited.gro -e etExcited.edr -g 
etExcited.log -pd &.


I've double checked my directory, all the file names are right and I've 
checked my nve.mdp file that I was running which I've posted below as well.


title   = production-dynamics
;warnings   = 10
cpp = /lib/cpp
;DEFINE = -DPOSRES
DEFINE = -DFLEXIBLE
;constraints = hbonds
;constraint_algorithm = shake
integrator = md
dt = 0.0005  ;
nsteps = 100  ;
;nstcomm   = 1 ; reset c.o.m motion
nstlist= 10
ns_type= grid
nstenergy  = 1  ; print energies
nstlog = 1  ; print to logfile
nstvout= 10   ; write velocities
nstxout= 10  ; collect data in fs (write coords)
;nstxtcout = 0   ; to print corrdinates to xtc trajectory
energygrps = Protein  HC4 r_42 r_46  r_50  r_52  r_68 r_69 r_70 SOL 
nstfout= 0

coulombtype= PME-Switch
fourierspacing = 0.12
pme_order  = 4
vdwtype= switch
rvdw   = 1.0
rlist  = 1.2
rcoulomb   = 1.0
pbc= xyz
;dispcorr   = Ener
continuation  = yes
;Berendsen temperature coupling is on
Tcoupl = no   ; temperature bath (yes, no)
;tau_t   = 0.1 
;tc-grps = system
;ref_t   = 300 
;Berendsen Pressure coupling is on

pcoupl  = no ; pressure bath (yes, no)
;pcoupltype  = isotropic
;tau_p   = 0.5
;compressibility = 4.5e-05
;ref_p   = 1.0
;Generate velocities  is on at 300
gen_vel = no   ; generate initial velocities
;gen_temp= 300.0 ; initial temperature
;gen_seed= 173529; random seeD
 
Lastly I used the command " grompp -f nve.mdp -c etExcited.gro -p 
etExcited.top -n res.ndx -o etExcited " to generate my .tpr file and I 
had no errors. To show that I have all the correct files in my directory 
 below I have printed the files from ls command.


nalysis   etSimulation.gro  #mdout.mdp.2#   posre_Protein.itp
E46QT50V.pdb   etSimulation.log  #mdout.mdp.3#   resE.ndx
em.mdp etSimulation.top  nohup.out   res.ndx
em.sh  etSimulation.trr  npt.mdp Save
etExcited.gro  et.tprnve.mdp state_prev.cpt
etExcited.top  machine   nvt.mdp stateSimulation.cpt
etExcited.tpr  mdout.mdp posre_Ion2.itp  topol.top
#etExcited.tpr.1#  #mdout.mdp.1# posre.itp

I have no idea why this stopped working for me, its just very recent.



If things used to work and now don't, and everything else is constant with 
respect to Gromacs (i.e. the version and specific installation you're using), 
then this is a problem for your sysadmin.  Sudden, unexpected changes mean 
someone is messing with your system.


-Justin


Thanks,

Taylor



*From:* Justin A. Lemkul 
*To:* Taylor Kaplan ; Discussion list for 
GROMACS users 

*Sent:* Friday, July 29, 2011 4:17 PM
*Subject:* Re: [gmx-users] cannot open file



Taylor Kaplan wrote:
 > Hi gromacs users,
 >
 >  I'm running a simulation but I keep getting the following errors.
 >
 > rogram mdrun_mpi, VERSION 4.5.1
 > Source code file: gmxfio.c, line: 537
 >
 > Can not open file:
 > taExcited.tpr
 > For more information and tips for troubleshooting, please check the 
GROMACS

 > website at http://www.gromacs.org/Documentation/Errors
 > ---
 >
 > "That Was Cool" (Beavis and Butthead)
 >
 > Error on node 0, will try to stop all the nodes
 > Halting parallel program mdrun_mpi on CPU 0 out of 4
 > 
--

 > MPI_ABORT was invoked on rank 0 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.
 > 
--

 >
 > gcq#182: "That Was Cool" (Beavis and Butthead)
 >
 > 
--

 > mpirun has exited due to process rank 0 with PID 9680 on
 > node compute-0-10 exiting without calling "finalize". This may
 > have caused other processes in the application to be
 > terminated by signals sent by mpirun (as reported here).
 >
 >  I don't know why gromacs cannot find my taExcited.tpr file, since it 
does exist in my working directory. However I noticed that when gromacs 
makes copies of my log files, its putting it in my home directory. T

Re: [gmx-users] bug in demux.pl

[Mark Abraham]

On 3/08/2011 8:19 AM, Swarnendu Tripathi wrote:

Hi,

I am not sure if this problem has been solved or discussed before.

I found that when I use the demux.pl  script it does 
not read the actual time from the .log file. This is for gromacs-4.0.7.


It reads the values fine, it just writes them in a limited-width format. 
You can find two occurrences of "%-8g" early in the script. Replace them 
with something like "%-20g" and things will work fine.


Thanks for the report.

Mark



These are few lines from one of the .log file for my replica exchnage 
simulation. For example in the second line below the actual time is  
"112.0" in ps. However, the demux.pl  script 
reads the time from the "Replica exchange" phrase as "1.1e+06". 
Same thing repeats in the following steps which is not right.

==
Step   Time Lambda
  28000 112.125000.0

   Energies (kJ/mol)
G96Bond   G96AngleProper Dih.  Improper Dih.  
LJ-14
1.58180e+031.49412e+037.67393e+023.69239e+02   
-2.78302e+02
 Coulomb-14LJ (SR)   Coulomb (SR)  Potential
Kinetic En.
   -2.39467e+008.01055e+01   -1.05400e+014.00143e+03
4.06353e+03

   Total EnergyTemperature Pressure (bar)
8.06496e+034.51898e+020.0e+00

Replica exchange at step 28000 time 1.1e+06
Repl 7 <-> 8  dE =  1.607e-01
Repl ex  01 x  23 x  45 x  67 x  89 x 10   11 x 
12   13 x 14   15 x 16   17 x 18   19 x 20   21 x 22   23
Repl pr.79   1.0   1.0   .85   1.0   
.98   .98   .90   1.0   1.0   .90


   Step   Time Lambda
  29000   116.125000.0

   Energies (kJ/mol)
G96Bond   G96AngleProper Dih.  Improper Dih.  
LJ-14
1.42588e+031.43257e+037.64199e+023.39804e+02   
-7.50646e+02
 Coulomb-14LJ (SR)   Coulomb (SR)  Potential
Kinetic En.
   -2.96172e+009.62908e+01   -1.95152e+013.28563e+03
4.16581e+03

   Total EnergyTemperature Pressure (bar)
7.45144e+034.63272e+020.0e+00

Replica exchange at step 29000 time 1.2e+06
Repl 8 <-> 9  dE =  2.023e-01
Repl ex  0 x  12 x  34 x  56 x  789   10 x 11   12 
x 13   14 x 15   16 x 17   18 x 19   20 x 21   22 x 23
Repl pr   1.0   .98   .66   1.0   .82   1.0   
.92   .93   .79   1.0   .96   1.0




As a result my "replica_index.xvg" looks like this which is not right.

1.19996e+06   119   22   12   10   215   23   177   19
831   150   136   20   1442   16   18
1.19996e+069   11   12   22   21   10   2357   178   
19130   156   13   14   2024   18   16
1.19997e+069   12   11   21   22   23   10758   17
1   19036   15   14   132   20   184   16
1.19997e+06   129   21   11   23   227   10851   
170   1963   14   152   13   20   18   164
1.19998e+06   12   219   23   11   2278   1051
0   176   19   1432   15   13   20   16   184
1.19998e+06   21   12   239   22   1187   1050
16   17   19   1423   13   15   16   204   18
1.19998e+06   21   23   12   2298   11   10750
61   19   172   14   133   16   154   20   18
1.1e+06   23   21   22   1298   10   11576
0   1912   17   13   14   1634   15   18   20
1.1e+06   23   22   219   12   1085   1167   
19021   13   17   16   1443   18   15   20
1.2e+0622   239   21   10   1258   116   19
720   131   16   174   14   183   20   15

=

Due to the problem in precision, time step repeats. The demux.pl 
 script should read the actual time from the previous 
step in the .log file.


-Swarnendu




-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] bug in demux.pl

[Swarnendu Tripathi]

I should have looked at the format when it writes the time but thanks for
your help.

-Swarnendu

On Tue, Aug 2, 2011 at 7:07 PM, Mark Abraham wrote:

>  On 3/08/2011 8:19 AM, Swarnendu Tripathi wrote:
>
> Hi,
>
> I am not sure if this problem has been solved or discussed before.
>
> I found that when I use the demux.pl script it does not read the actual
> time from the .log file. This is for gromacs-4.0.7.
>
>
> It reads the values fine, it just writes them in a limited-width format.
> You can find two occurrences of "%-8g" early in the script. Replace them
> with something like "%-20g" and things will work fine.
>
> Thanks for the report.
>
> Mark
>
>
>
> These are few lines from one of the .log file for my replica exchnage
> simulation. For example in the second line below the actual time is
> "112.0" in ps. However, the demux.pl script reads the time from the
> "Replica exchange" phrase as "1.1e+06". Same thing repeats in the
> following steps which is not right.
> ==
> Step   Time Lambda
>   28000   112.125000.0
>
>Energies (kJ/mol)
> G96Bond   G96AngleProper Dih.  Improper Dih.  LJ-14
> 1.58180e+031.49412e+037.67393e+023.69239e+02   -2.78302e+02
>  Coulomb-14LJ (SR)   Coulomb (SR)  PotentialKinetic En.
>-2.39467e+008.01055e+01   -1.05400e+014.00143e+034.06353e+03
>Total EnergyTemperature Pressure (bar)
> 8.06496e+034.51898e+020.0e+00
>
> Replica exchange at step 28000 time 1.1e+06
> Repl 7 <-> 8  dE =  1.607e-01
> Repl ex  01 x  23 x  45 x  67 x  89 x 10   11 x 12   13
> x 14   15 x 16   17 x 18   19 x 20   21 x 22   23
> Repl pr.79   1.0   1.0   .85   1.0   .98
> .98   .90   1.0   1.0   .90
>
>Step   Time Lambda
>   29000   116.125000.0
>
>Energies (kJ/mol)
> G96Bond   G96AngleProper Dih.  Improper Dih.  LJ-14
> 1.42588e+031.43257e+037.64199e+023.39804e+02   -7.50646e+02
>  Coulomb-14LJ (SR)   Coulomb (SR)  PotentialKinetic En.
>-2.96172e+009.62908e+01   -1.95152e+013.28563e+034.16581e+03
>Total EnergyTemperature Pressure (bar)
> 7.45144e+034.63272e+020.0e+00
>
> Replica exchange at step 29000 time 1.2e+06
> Repl 8 <-> 9  dE =  2.023e-01
> Repl ex  0 x  12 x  34 x  56 x  789   10 x 11   12 x
> 13   14 x 15   16 x 17   18 x 19   20 x 21   22 x 23
> Repl pr   1.0   .98   .66   1.0   .82   1.0
> .92   .93   .79   1.0   .96   1.0
>
>
> 
>
> As a result my "replica_index.xvg" looks like this which is not right.
>
> 
> 1.19996e+06   119   22   12   10   215   23   177   198
> 31   150   136   20   1442   16   18
> 1.19996e+069   11   12   22   21   10   2357   178   19
> 130   156   13   14   2024   18   16
> 1.19997e+069   12   11   21   22   23   10758   171
> 19036   15   14   132   20   184   16
> 1.19997e+06   129   21   11   23   227   10851   17
> 0   1963   14   152   13   20   18   164
> 1.19998e+06   12   219   23   11   2278   10510
> 176   19   1432   15   13   20   16   184
> 1.19998e+06   21   12   239   22   1187   10501
> 6   17   19   1423   13   15   16   204   18
> 1.19998e+06   21   23   12   2298   11   107506
> 1   19   172   14   133   16   154   20   18
> 1.1e+06   23   21   22   1298   10   115760
> 1912   17   13   14   1634   15   18   20
> 1.1e+06   23   22   219   12   1085   1167   19
> 021   13   17   16   1443   18   15   20
> 1.2e+0622   239   21   10   1258   116   197
> 20   131   16   174   14   183   20   15
>
> =
>
> Due to the problem in precision, time step repeats. The demux.pl script
> should read the actual time from the previous step in the .log file.
>
> -Swarnendu
>
>
>
>
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> Can't post? 

Re: [gmx-users] spatial distribution function in a binary solvent mixture

[dbiswal]

 Thanks Chris. I presume g_sdf won't be helpful for my system.

> I think that you have a misconception about what g_spatial does. For a
> system with many type A and many type B, you need to average over all
> of one type as the central solute to compute an rdf, and perhaps that
> is what you want for your sdf. g_spatial, however, does not do any
> fitting. g_sdf did that. You can obtain an old version of the source
> (4.0.7 should have it I think) if you want to use g_sdf. I have never
> used g_sdf myself.
>
> In case that doesn't answer your question, then let me make one more
> point:
>
> g_spatial requires that a bin exists for every count. Thus either (a)
> find a large memory node somewhere or (b) pre-process your trajectory
> using trjorder to order the solvent molecules and then keep only the N
> closest, then run g_spatial. But again, I think that this will not
> provide what you are looking for but is likely to give you a smear
> over your box.
>
> Chris.
>
> -- original message --
>
>Dear all,
>
>   I'm working with a binary solvent mixture containing 2000 molecules
> (1800
>   type-A + 200 type-B). Both the types of solvent molecules have similar
> structure (they are both diatomic molecules) except the polarity.I'm
> trying to calculate sdf of type-B solvent molecules. I followed the
> step-by-step instructions from the manual using g_spatial.
>   I'm trying to reduce the bin width to 0.05A. The *.cube file generated
> with a bin width of 0.09A is already 4.9GB in size. As I'm more
> interested in the first solvation shell around the type-B solvent
> molecule, I was wondering if I could find a way to control the maximum
> radius of the sphere around central molecule (center of coordinates?)?
> Also, can anyone please let me know how g_spatial deals with the angular
> part of sdf?
>
> (I've searched a lot but could not gather enough meaningful information)
>
> regards,
> Debasmita
>
>
>
>
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
>
>


-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

[gmx-users] Problem with temperature and stochastic dynamics in FEP

[Sanku M]

Hi,
  I am trying to calculate the solvation free energy using thermodynamic 
integration(TI) method . I am using gromacs 4.0.7 . But, I am having a problem 
in getting accurate average temperature .
The following is what I did:
   I found that when doing TI, grompp recommends using 'sd' ( stochastic 
dynamics ) as integrator ( provides an warning that for decoupled system, it is 
better to use sd in stead of md ).
Also I went through Justin Lemkul's tutorial which also recommends using 'sd' 
for TI method.   So, I tried to use sd integrator. 

I wanted to run NVT simulation at 300 K.
I also found from manual and also from Justin's website on FEP, sd integrator 
implicitly controls temperature and so there is no need to specify a 
thermostat. 
So, I did not specify any thermostat ( I wrote tcoupl = No ). But, 
unfortunately, I found that after a long simulation  , the average temperature 
actually goes to 303 K in stead of 300 K( the desired temperature). This 
happened for all Lambda values where the average temperature turned out to be 
303 K. The .mdp file is pasted at the end of the email and I felt I am using 
reasonable cutoffs and PME as electrostatics.

Here is the output from g_energy command to get the temperature. 
 
 Statistics over 401 steps [ 0. thru 8000.0005 ps ], 1 data sets
All averages are exact over 401 steps

Energy  Average   RMSD Fluct.  Drift  Tot-Drift
---
Temperature 303.3334.780844.78082 5.46937e-06   0.043755
Heat Capacity Cv:  12.4764 J/mol K (factor = 0.00024841)


I tried three  further tests:
a) I removed tcoupl = No option . But, it is giving same 303K as average 
temperature when using sd.

a)  I tried to specify the Nose-Hoover thermostat along with sd. But still, I 
found that when using sd, the average  temperature is going to 303 K in stead 
of 
300 K.

c)  Finally, I tried to overlook grompp warning and went ahead and used 'md' 
along with Nose-Hoover thermostat. This time I found the right average 
temperature 300 K is being achieved.B
But, I understand that for a decoupled system like here, I need to use a method 
like stochastic dynamics. 

 But , I was wondering why it is reaching 303 K in stead of 300 K when using sd 
 but 300K is achieved when using md. I looked at the mailing list where people 
had issues with sd and temperature control. But, I could not find a good 
solution. So, any help on the right protocol will be really appreciated.


Here is my .mdp file.

; RUN CONTROL PARAMETERS
integrator   = sd
; Start time and timestep in ps
tinit= 0.0
dt   = 0.002
nsteps   = 400
; For exact run continuation or redoing part of a run
; Part index is updated automatically on checkpointing (keeps files separate)
simulation_part  = 1
init_step= 0
; mode for center of mass motion removal
comm-mode= Linear
; number of steps for center of mass motion removal
nstcomm  = 1
; group(s) for center of mass motion removal
comm-grps=
; LANGEVIN DYNAMICS OPTIONS
; Friction coefficient (amu/ps) and random seed
bd-fric  = 0
ld-seed  = 1993



; OUTPUT CONTROL OPTIONS
; OUTPUT CONTROL OPTIONS
; Output frequency for coords (x), velocities (v) and forces (f)
nstxout  = 1
nstvout  = 1
nstfout  = 0
; Output frequency for energies to log file and energy file
nstlog   = 1000
nstenergy= 10
; Output frequency and precision for xtc file
nstxtcout= 1000
xtc-precision= 1000
; This selects the subset of atoms for the xtc file. You can
; select multiple groups. By default all atoms will be written.
xtc-grps =
; Selection of energy groups
energygrps   =

; NEIGHBORSEARCHING PARAMETERS
; nblist update frequency
nstlist  = 10
; ns algorithm (simple or grid)
ns_type  = grid
; Periodic boundary conditions: xyz, no, xy
pbc  = xyz
periodic_molecules   = no
; nblist cut-off
rlist= 1.4

; OPTIONS FOR ELECTROSTATICS AND VDW
; Method for doing electrostatics
coulombtype  = pme
rcoulomb-switch  = 0
rcoulomb = 1.4
; Relative dielectric constant for the medium and the reaction field
epsilon_r= 1
epsilon_rf   = 1
; Method for doing Van der Waals
epsilon_rf   = 1
; Method for doing Van der Waals
vdw-type = shift
; cut-off lengths
rvdw-switch  = 1.0
rvdw = 1.4
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = EnerPres
; Extension of the potential lookup tables beyond the cut-off
table-extension  = 1
; Seperat

Re: [gmx-users] Problem with temperature and stochastic dynamics in FEP

[Justin A. Lemkul]

Sanku M wrote:

Hi,
  I am trying to calculate the solvation free energy using thermodynamic 
integration(TI) method . I am using gromacs 4.0.7 . But, I am having a 
problem in getting accurate average temperature .

The following is what I did:
   I found that when doing TI, grompp recommends using 'sd' ( stochastic 
dynamics ) as integrator ( provides an warning that for decoupled 
system, it is better to use sd in stead of md ).
Also I went through Justin Lemkul's tutorial which also recommends using 
'sd' for TI method.   So, I tried to use sd integrator. 



The tutorial uses BAR for calculation free energy, not TI.  Though the protocols 
are similar, there is a difference, FYI.



I wanted to run NVT simulation at 300 K.
I also found from manual and also from Justin's website on FEP, sd 
integrator implicitly controls temperature and so there is no need to 
specify a thermostat. So, I did not specify any thermostat ( I wrote 
tcoupl = No ). But, unfortunately, I found that after a long simulation 
 , the average temperature actually goes to 303 K in stead of 300 K( the 
desired temperature). This happened for all Lambda values where the 
average temperature turned out to be 303 K. The .mdp file is pasted at 
the end of the email and I felt I am using reasonable cutoffs and PME as 
electrostatics.


Here is the output from g_energy command to get the temperature. 
 
 Statistics over 401 steps [ 0. thru 8000.0005 ps ], 1 data sets

All averages are exact over 401 steps

Energy  Average   RMSD Fluct.  Drift 
 Tot-Drift

---
Temperature 303.3334.780844.78082 5.46937e-06   
0.043755

Heat Capacity Cv:  12.4764 J/mol K (factor = 0.00024841)


I tried three  further tests:
a) I removed tcoupl = No option . But, it is giving same 303K as average 
temperature when using sd.


a)  I tried to specify the Nose-Hoover thermostat along with sd. But 
still, I found that when using sd, the average  temperature is going to 
303 K in stead of 300 K.




From the manual description of the sd integrator: "The parameter tcoupl is 
ignored."  This explains (a) and (b).


c)  Finally, I tried to overlook grompp warning and went ahead and used 
'md' along with Nose-Hoover thermostat. This time I found the right 
average temperature 300 K is being achieved.B
But, I understand that for a decoupled system like here, I need to use a 
method like stochastic dynamics. 

 But , I was wondering why it is reaching 303 K in stead of 300 K when 
using sd  but 300K is achieved when using md. I looked at the mailing 
list where people had issues with sd and temperature control. But, I 
could not find a good solution. So, any help on the right protocol will 
be really appreciated.




What type of equilibration did you do prior to the data collection?  If your 
system isn't sampling the desired ensemble, then you shouldn't proceed.


-Justin



Here is my .mdp file.

; RUN CONTROL PARAMETERS
integrator   = sd
; Start time and timestep in ps
tinit= 0.0
dt   = 0.002
nsteps   = 400
; For exact run continuation or redoing part of a run
; Part index is updated automatically on checkpointing (keeps files 
separate)

simulation_part  = 1
init_step= 0
; mode for center of mass motion removal
comm-mode= Linear
; number of steps for center of mass motion removal
nstcomm  = 1
; group(s) for center of mass motion removal
comm-grps=
; LANGEVIN DYNAMICS OPTIONS
; Friction coefficient (amu/ps) and random seed
bd-fric  = 0
ld-seed  = 1993



; OUTPUT CONTROL OPTIONS
; OUTPUT CONTROL OPTIONS
; Output frequency for coords (x), velocities (v) and forces (f)
nstxout  = 1
nstvout  = 1
nstfout  = 0
; Output frequency for energies to log file and energy file
nstlog   = 1000
nstenergy= 10
; Output frequency and precision for xtc file
nstxtcout= 1000
xtc-precision= 1000
; This selects the subset of atoms for the xtc file. You can
; select multiple groups. By default all atoms will be written.
xtc-grps =
; Selection of energy groups
energygrps   =

; NEIGHBORSEARCHING PARAMETERS
; nblist update frequency
nstlist  = 10
; ns algorithm (simple or grid)
ns_type  = grid
; Periodic boundary conditions: xyz, no, xy
pbc  = xyz
periodic_molecules   = no
; nblist cut-off
rlist= 1.4

; OPTIONS FOR ELECTROSTATICS AND VDW
; Method for doing electrostatics
coulombtype  = pme
rcoulomb-switch  = 0
rcoulomb = 1.4
; Relative dielectric constant for the medium and the reaction field
epsilon_r= 

Re: [gmx-users] Problem with temperature and stochastic dynamics in FEP

[Sanku M]

Hi Justin,
  I first performed two minimization using steep and l-bfgs as the method on 
the 
solvated system for each lambda. Then I ran a 500ps NVT simulation on minimized 
system ( for each lambda) where essentially same .mdp file was used except the 
following change: 

   gen-vel  = yes gen-temp= 280
 gen-seed= -1


So, essentially, I heated the minimized system up from 280 K with ref_t 300K 
using sd integrator. But, here also I found the average temperature goes to 303 
K( instead of 300 K) .
One thing I don't understand that why using md as integrator gives the desired 
average temperature( when using Nose-Hoover thermostat) while the sd integrator 
does not.
Any further suggestion will be appreciated.

Jagannath



From: Justin A. Lemkul 
To: Discussion list for GROMACS users 
Sent: Tue, August 2, 2011 8:08:20 PM
Subject: Re: [gmx-users] Problem with temperature and stochastic dynamics in FEP



Sanku M wrote:
> Hi,
>   I am trying to calculate the solvation free energy using thermodynamic 
>integration(TI) method . I am using gromacs 4.0.7 . But, I am having a problem 
>in getting accurate average temperature .
> The following is what I did:
>I found that when doing TI, grompp recommends using 'sd' ( stochastic 
>dynamics ) as integrator ( provides an warning that for decoupled system, it 
>is 
>better to use sd in stead of md ).
> Also I went through Justin Lemkul's tutorial which also recommends using 'sd' 
>for TI method.   So, I tried to use sd integrator. 
>

The tutorial uses BAR for calculation free energy, not TI.  Though the 
protocols 
are similar, there is a difference, FYI.

> I wanted to run NVT simulation at 300 K.
> I also found from manual and also from Justin's website on FEP, sd integrator 
>implicitly controls temperature and so there is no need to specify a 
>thermostat. 
>So, I did not specify any thermostat ( I wrote tcoupl = No ). But, 
>unfortunately, I found that after a long simulation  , the average temperature 
>actually goes to 303 K in stead of 300 K( the desired temperature). This 
>happened for all Lambda values where the average temperature turned out to be 
>303 K. The .mdp file is pasted at the end of the email and I felt I am using 
>reasonable cutoffs and PME as electrostatics.
> 
> Here is the output from g_energy command to get the temperature.   Statistics 
>over 401 steps [ 0. thru 8000.0005 ps ], 1 data sets
> All averages are exact over 401 steps
> 
> Energy  Average   RMSD Fluct.  Drift  
Tot-Drift
> 
---
> Temperature 303.3334.780844.78082 5.46937e-06   
>0.043755
> Heat Capacity Cv:  12.4764 J/mol K (factor = 0.00024841)
> 
> 
> I tried three  further tests:
> a) I removed tcoupl = No option . But, it is giving same 303K as average 
>temperature when using sd.
> 
> a)  I tried to specify the Nose-Hoover thermostat along with sd. But still, I 
>found that when using sd, the average  temperature is going to 303 K in stead 
>of 
>300 K.
> 

>From the manual description of the sd integrator: "The parameter tcoupl is 
ignored."  This explains (a) and (b).

> c)  Finally, I tried to overlook grompp warning and went ahead and used 'md' 
>along with Nose-Hoover thermostat. This time I found the right average 
>temperature 300 K is being achieved.B
> But, I understand that for a decoupled system like here, I need to use a 
> method 
>like stochastic dynamics. 
>
>  But , I was wondering why it is reaching 303 K in stead of 300 K when using 
>sd  but 300K is achieved when using md. I looked at the mailing list where 
>people had issues with sd and temperature control. But, I could not find a 
>good 
>solution. So, any help on the right protocol will be really appreciated.
> 

What type of equilibration did you do prior to the data collection?  If your 
system isn't sampling the desired ensemble, then you shouldn't proceed.

-Justin

> 
> Here is my .mdp file.
> 
> ; RUN CONTROL PARAMETERS
> integrator   = sd
> ; Start time and timestep in ps
> tinit= 0.0
> dt   = 0.002
> nsteps   = 400
> ; For exact run continuation or redoing part of a run
> ; Part index is updated automatically on checkpointing (keeps files separate)
> simulation_part  = 1
> init_step= 0
> ; mode for center of mass motion removal
> comm-mode= Linear
> ; number of steps for center of mass motion removal
> nstcomm  = 1
> ; group(s) for center of mass motion removal
> comm-grps=
> ; LANGEVIN DYNAMICS OPTIONS
> ; Friction coefficient (amu/ps) and random seed
> bd-fric  = 0
> ld-seed  = 1993
> 
> 
> 
> ; OUTPUT CONTROL OPTIONS
> ; OUTPUT CONTROL OPTIONS
> ; Output frequency for coords (x), velocit

Re: [gmx-users] Problem with temperature and stochastic dynamics in FEP

[Justin A. Lemkul]

Sanku M wrote:

Hi Justin,
  I first performed two minimization using steep and l-bfgs as the 
method on the solvated system for each lambda. Then I ran a 500ps NVT 
simulation on minimized system ( for each lambda) where essentially same 
.mdp file was used except the following change: 


   gen-vel  = yes
 gen-temp= 280
 gen-seed= -1

So, essentially, I heated the minimized system up from 280 K with ref_t 
300K using sd integrator. But, here also I found the average temperature 
goes to 303 K( instead of 300 K) .


Unless you did simulated annealing in between, you didn't heat the system, you 
suddenly jolted the temperature, which is probably not stable.  Maybe the 
integrator/Langevin thermostat does not react well to these conditions. 
Equilibrate under the conditions you wish to collect data.  Your data collection 
uses NPT, so running NVT at a different temperature and suddenly changing the 
temperature and introducing pressure coupling is not an appropriate procedure. 
If nothing else, your data will be unreliable, even if your temperature was correct.


One thing I don't understand that why using md as integrator gives the 
desired average temperature( when using Nose-Hoover thermostat) while 
the sd integrator does not.


I can only hazard a guess that for some reason Nose-Hoover can arrive at the 
correct temperature more quickly, but that seems uncharacteristic for it given 
the nature of its fluctuations.  Perhaps someone with more knowledge of the 
algorithms and code can comment.  But until you've established a properly 
equilibrated system under appropriate conditions, it's not worth investigating 
these differences much further.


-Justin


Any further suggestion will be appreciated.

Jagannath

*From:* Justin A. Lemkul 
*To:* Discussion list for GROMACS users 
*Sent:* Tue, August 2, 2011 8:08:20 PM
*Subject:* Re: [gmx-users] Problem with temperature and stochastic 
dynamics in FEP




Sanku M wrote:
 > Hi,
 >  I am trying to calculate the solvation free energy using 
thermodynamic integration(TI) method . I am using gromacs 4.0.7 . But, I 
am having a problem in getting accurate average temperature .

 > The following is what I did:
 >I found that when doing TI, grompp recommends using 'sd' ( 
stochastic dynamics ) as integrator ( provides an warning that for 
decoupled system, it is better to use sd in stead of md ).
 > Also I went through Justin Lemkul's tutorial which also recommends 
using 'sd' for TI method.  So, I tried to use sd integrator.


The tutorial uses BAR for calculation free energy, not TI.  Though the 
protocols are similar, there is a difference, FYI.


 > I wanted to run NVT simulation at 300 K.
 > I also found from manual and also from Justin's website on FEP, sd 
integrator implicitly controls temperature and so there is no need to 
specify a thermostat. So, I did not specify any thermostat ( I wrote 
tcoupl = No ). But, unfortunately, I found that after a long simulation  
, the average temperature actually goes to 303 K in stead of 300 K( the 
desired temperature). This happened for all Lambda values where the 
average temperature turned out to be 303 K. The .mdp file is pasted at 
the end of the email and I felt I am using reasonable cutoffs and PME as 
electrostatics.

 >
 > Here is the output from g_energy command to get the temperature.  
Statistics over 401 steps [ 0. thru 8000.0005 ps ], 1 data sets

 > All averages are exact over 401 steps
 >
 > Energy  Average  RMSDFluct.  Drift  
Tot-Drift
 > 
---
 > Temperature303.3334.780844.78082 5.46937e-06  
0.043755

 > Heat Capacity Cv:  12.4764 J/mol K (factor = 0.00024841)
 >
 >
 > I tried three  further tests:
 > a) I removed tcoupl = No option . But, it is giving same 303K as 
average temperature when using sd.

 >
 > a)  I tried to specify the Nose-Hoover thermostat along with sd. But 
still, I found that when using sd, the average  temperature is going to 
303 K in stead of 300 K.

 >

 From the manual description of the sd integrator: "The parameter tcoupl 
is ignored."  This explains (a) and (b).


 > c)  Finally, I tried to overlook grompp warning and went ahead and 
used 'md' along with Nose-Hoover thermostat. This time I found the right 
average temperature 300 K is being achieved.B
 > But, I understand that for a decoupled system like here, I need to 
use a method like stochastic dynamics.
 >  But , I was wondering why it is reaching 303 K in stead of 300 K 
when using sd  but 300K is achieved when using md. I looked at the 
mailing list where people had issues with sd and temperature control. 
But, I could not find a good solution. So, any help on the right 
protocol will be really appreciated.

 >

What type of equilibration did you

Re: [gmx-users] Problem with temperature and stochastic dynamics in FEP

[Justin A. Lemkul]

Justin A. Lemkul wrote:



Sanku M wrote:

Hi Justin,
  I first performed two minimization using steep and l-bfgs as the 
method on the solvated system for each lambda. Then I ran a 500ps NVT 
simulation on minimized system ( for each lambda) where essentially 
same .mdp file was used except the following change:

   gen-vel  = yes
 gen-temp= 280
 gen-seed= -1

So, essentially, I heated the minimized system up from 280 K with 
ref_t 300K using sd integrator. But, here also I found the average 
temperature goes to 303 K( instead of 300 K) .


Unless you did simulated annealing in between, you didn't heat the 
system, you suddenly jolted the temperature, which is probably not 
stable.  Maybe the integrator/Langevin thermostat does not react well to 
these conditions. Equilibrate under the conditions you wish to collect 
data.  Your data collection uses NPT, so running NVT at a different 
temperature and suddenly changing the temperature and introducing 
pressure coupling is not an appropriate procedure. If nothing else, your 
data will be unreliable, even if your temperature was correct.




Sorry, misread the previous .mdp settings.  You're still using NVT, but the 
point about the temperature is still valid here.


-Justin

--


Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.

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

[gmx-users] spatial distribution function in a binary solvent mixture

[chris . neale]

Here's the totally wasteful way for you to get what you want:

1. for each molecule in type A: trjorder everything around that  
molecule and output only the closest N atoms of type A and the closest  
M atoms of type B. Ensure that the coordinate order is: your central  
molecule of type A is first, the remaining type A second, and the type  
B third.


2. concatenate all of these trajectories into one single large trajectory

3. trjconv -fit trans+rot , fitting only to the single central molecule.

4a. g_spatial for the non-central type A
4b. g_spatial for type B

* now repeat switching types A and B in the above steps to get the  
SDFs arouns type B.


A much better idea is for you to modify g_spatial to do this all  
within one loop. But if you don't know C then you are out of luck there.


g_sdf may not work because you are underdetermined (there is a  
symmetry about the long axis that you can not specify). But perhaps  
you could use that to your advantage and select some totally unrelated  
atom as the third atom for the fit. My best guess is that this would  
be no worse than g_spatial.


A final option is for you to use g_mindist to output all of the  
contacts and then run over that with your own script to process the  
data into an SDF and output it in cube format.


Chris.

-- original message --

 Thanks Chris. I presume g_sdf won't be helpful for my system.

[Hide Quoted Text]
I think that you have a misconception about what g_spatial does. For a
system with many type A and many type B, you need to average over all
of one type as the central solute to compute an rdf, and perhaps that
is what you want for your sdf. g_spatial, however, does not do any
fitting. g_sdf did that. You can obtain an old version of the source
(4.0.7 should have it I think) if you want to use g_sdf. I have never
used g_sdf myself.

In case that doesn't answer your question, then let me make one more
point:

g_spatial requires that a bin exists for every count. Thus either (a)
find a large memory node somewhere or (b) pre-process your trajectory
using trjorder to order the solvent molecules and then keep only the N
closest, then run g_spatial. But again, I think that this will not
provide what you are looking for but is likely to give you a smear
over your box.

Chris.

-- original message --

Dear all,

   I'm working with a binary solvent mixture containing 2000 molecules
(1800
   type-A + 200 type-B). Both the types of solvent molecules have similar
structure (they are both diatomic molecules) except the polarity.I'm
trying to calculate sdf of type-B solvent molecules. I followed the
step-by-step instructions from the manual using g_spatial.
   I'm trying to reduce the bin width to 0.05A. The *.cube file generated
with a bin width of 0.09A is already 4.9GB in size. As I'm more
interested in the first solvation shell around the type-B solvent
molecule, I was wondering if I could find a way to control the maximum
radius of the sphere around central molecule (center of coordinates?)?
Also, can anyone please let me know how g_spatial deals with the angular
part of sdf?

(I've searched a lot but could not gather enough meaningful information)

regards,
Debasmita

--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Problem with temperature and stochastic dynamics in FEP

[Sanku M]

Hi Justin,
  So, do you suggest that after minimization, I should generate the velocity at 
300K instead ?
i.,e for equilibration,should following the set-up?

gen-vel  = yes
>>  gen-temp= 300
>>  gen-seed= -1




From: Justin A. Lemkul 
To: Discussion list for GROMACS users 
Sent: Tue, August 2, 2011 8:26:16 PM
Subject: Re: [gmx-users] Problem with temperature and stochastic dynamics in FEP



Justin A. Lemkul wrote:
> 
> 
> Sanku M wrote:
>> Hi Justin,
>>   I first performed two minimization using steep and l-bfgs as the method on 
>>the solvated system for each lambda. Then I ran a 500ps NVT simulation on 
>>minimized system ( for each lambda) where essentially same .mdp file was used 
>>except the following change:
>>gen-vel  = yes
>>  gen-temp= 280
>>  gen-seed= -1
>> 
>> So, essentially, I heated the minimized system up from 280 K with ref_t 300K 
>>using sd integrator. But, here also I found the average temperature goes to 
>>303 
>>K( instead of 300 K) .
> 
> Unless you did simulated annealing in between, you didn't heat the system, 
> you 
>suddenly jolted the temperature, which is probably not stable.  Maybe the 
>integrator/Langevin thermostat does not react well to these conditions. 
>Equilibrate under the conditions you wish to collect data.  Your data 
>collection 
>uses NPT, so running NVT at a different temperature and suddenly changing the 
>temperature and introducing pressure coupling is not an appropriate procedure. 
>If nothing else, your data will be unreliable, even if your temperature was 
>correct.
> 

Sorry, misread the previous .mdp settings.  You're still using NVT, but the 
point about the temperature is still valid here.

-Justin

-- 

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin


-- gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search 
before posting!
Please don't post (un)subscribe requests to the list. Use the www interface or 
send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Problem with temperature and stochastic dynamics in FEP

[Justin A. Lemkul]

Sanku M wrote:

Hi Justin,
  So, do you suggest that after minimization, I should generate the 
velocity at 300K instead ?

i.,e for equilibration,should following the set-up?



You should always equilibrate under the desired conditions.  Never make sudden 
changes, or else you're going to damage the stability of the system.


-Justin


gen-vel  = yes

>  gen-temp= 300
>  gen-seed= -1



*From:* Justin A. Lemkul 
*To:* Discussion list for GROMACS users 
*Sent:* Tue, August 2, 2011 8:26:16 PM
*Subject:* Re: [gmx-users] Problem with temperature and stochastic 
dynamics in FEP




Justin A. Lemkul wrote:
 >
 >
 > Sanku M wrote:
 >> Hi Justin,
 >>  I first performed two minimization using steep and l-bfgs as the 
method on the solvated system for each lambda. Then I ran a 500ps NVT 
simulation on minimized system ( for each lambda) where essentially same 
.mdp file was used except the following change:

 >>gen-vel  = yes
 >>  gen-temp= 280
 >>  gen-seed= -1
 >>
 >> So, essentially, I heated the minimized system up from 280 K with 
ref_t 300K using sd integrator. But, here also I found the average 
temperature goes to 303 K( instead of 300 K) .

 >
 > Unless you did simulated annealing in between, you didn't heat the 
system, you suddenly jolted the temperature, which is probably not 
stable.  Maybe the integrator/Langevin thermostat does not react well to 
these conditions. Equilibrate under the conditions you wish to collect 
data.  Your data collection uses NPT, so running NVT at a different 
temperature and suddenly changing the temperature and introducing 
pressure coupling is not an appropriate procedure. If nothing else, your 
data will be unreliable, even if your temperature was correct.

 >

Sorry, misread the previous .mdp settings.  You're still using NVT, but 
the point about the temperature is still valid here.


-Justin

-- 

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin


-- gmx-users mailing listgmx-users@gromacs.org 


http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the www 
interface or send it to gmx-users-requ...@gromacs.org 
.

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


--


Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.

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

Re: [gmx-users] Problem with temperature and stochastic dynamics in FEP

[Sanku M]

OKK, But, before, I have tried using md as integrator where I start the initial 
velocity as 280K ( which makes sense as the minimized configuration is 
essentially at low temperature) and have used Nose-hoover at 300K and I never 
had problem in getting the desired 300K after running a short NVT simulation as 
equilabration.
The problem arises when sd integrator is used.
In fact, I  just tried to generate the initial velocity at 300K along with sd 
integrator. I got the 303 K as average temperature again. The same simulation 
using md integrator always gives rise to 300K as average temperature( when used 
in conjunction with Nose-Hoover ) , even if I generate velocity at 280 K or 
300K 
.

I guess that the initial velocity-generation temperature should not be an issue 
here. The sd integrator in contrast to md integrator may be the problem.






From: Justin A. Lemkul 
To: Discussion list for GROMACS users 
Sent: Tue, August 2, 2011 8:40:21 PM
Subject: Re: [gmx-users] Problem with temperature and stochastic dynamics in FEP



Sanku M wrote:
> Hi Justin,
>   So, do you suggest that after minimization, I should generate the velocity 
> at 
>300K instead ?
> i.,e for equilibration,should following the set-up?
> 

You should always equilibrate under the desired conditions.  Never make sudden 
changes, or else you're going to damage the stability of the system.

-Justin

> gen-vel  = yes
>> >  gen-temp= 300
>> >  gen-seed= -1
> 
> 
> *From:* Justin A. Lemkul 
> *To:* Discussion list for GROMACS users 
> *Sent:* Tue, August 2, 2011 8:26:16 PM
> *Subject:* Re: [gmx-users] Problem with temperature and stochastic dynamics 
> in 
>FEP
> 
> 
> 
> Justin A. Lemkul wrote:
>  >
>  >
>  > Sanku M wrote:
>  >> Hi Justin,
>  >>  I first performed two minimization using steep and l-bfgs as the method 
> on 
>the solvated system for each lambda. Then I ran a 500ps NVT simulation on 
>minimized system ( for each lambda) where essentially same .mdp file was used 
>except the following change:
>  >>gen-vel  = yes
>  >>  gen-temp= 280
>  >>  gen-seed = -1
>  >>
>  >> So, essentially, I heated the minimized system up from 280 K with ref_t 
>300K using sd integrator. But, here also I found the average temperature goes 
>to 
>303 K( instead of 300 K) .
>  >
>  > Unless you did simulated annealing in between, you didn't heat the system, 
>you suddenly jolted the temperature, which is probably not stable.  Maybe the 
>integrator/Langevin thermostat does not react well to these conditions. 
>Equilibrate under the conditions you wish to collect data.  Your data 
>collection 
>uses NPT, so running NVT at a different temperature and suddenly changing the 
>temperature and introducing pressure coupling is not an appropriate procedure. 
>If nothing else, your data will be unreliable, even if your temperature was 
>correct.
>  >
> 
> Sorry, misread the previous .mdp settings.  You're still using  NVT, but the 
>point about the temperature is still valid here.
> 
> -Justin
> 
> -- 
> 
> Justin A. Lemkul
> Ph.D. Candidate
> ICTAS Doctoral Scholar
> MILES-IGERT Trainee
> Department of Biochemistry
> Virginia Tech
> Blacksburg, VA
> jalemkul[at]vt.edu | (540) 231-9080
> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
> 
> 
> -- gmx-users mailing listgmx-users@gromacs.org 
>
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at 
>http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> Please don't post (un)subscribe requests to the list. Use the www interface 
> or 
>send it to gmx-users-requ...@gromacs.org 
>.
> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

-- 

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia  Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin


-- gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search 
before posting!
Please don't post (un)subscribe requests to the list. Use the www interface or 
send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Blacksburg, VA
> jalemkul[at]vt.edu | (540) 231-9080
> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
> 
> 
> -- gmx-users mailing listgmx-users@gromacs.org 
>
> http://lists.gromacs.org/mailman/listinf

Re: [gmx-users] Problem with temperature and stochastic dynamics in FEP

[Justin A. Lemkul]

Sanku M wrote:
OKK, But, before, I have tried using md as integrator where I start the 
initial velocity as 280K ( which makes sense as the minimized 
configuration is essentially at low temperature) and have used 
Nose-hoover at 300K and I never had problem in getting the desired 300K 
after running a short NVT simulation as equilabration.

The problem arises when sd integrator is used.
In fact, I  just tried to generate the initial velocity at 300K along 
with sd integrator. I got the 303 K as average temperature again. The 
same simulation using md integrator always gives rise to 300K as average 
temperature( when used in conjunction with Nose-Hoover ) , even if I 
generate velocity at 280 K or 300K .


I guess that the initial velocity-generation temperature should not be 
an issue here. The sd integrator in contrast to md integrator may be the 
problem.





I have never had such a problem when using sd.  I can always achieve the target 
temperature within 0.5 K, if not exactly, so I do not believe there is an 
inherent problem with sd.  Your first post indicated that your statistics were 
collected over the entire timeframe, which will not be correct for reasons I 
have already described.  The initial frames will represent unequilibrated data 
that may skew the averages.  What if you analyze only the last half of the time? 
 The last 1 ns?


-Justin

--


Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.

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

Re: [gmx-users] Problem with temperature and stochastic dynamics in FEP

[Mark Abraham]

On 03/08/11, Sanku M   wrote:
> 
> 
> 
> 
> 
> 
> OKK, But, before, I have tried using md as integrator where I start the 
> initial velocity as 280K ( which makes sense as the minimized configuration 
> is essentially at low temperature) and have used Nose-hoover at 300K and I 
> never had problem in getting the desired 300K after running a short NVT 
> simulation as equilabration.
> The problem arises when sd integrator is used.
> In fact, I  just tried to generate the initial velocity at 300K along with sd 
> integrator. I got the 303 K as average temperature again. The same simulation 
> using md integrator always gives rise to 300K as average temperature( when 
> used in conjunction with Nose-Hoover ) , even if I generate velocity at 280 K 
> or 300K .
> 
> 
> I guess that the initial velocity-generation temperature should not be an 
> issue here. The sd integrator in contrast to md integrator may be the
>  problem.
> 
> 
> 
> 
> 
> 

Or your use of it... check out manual 7.3.3 and how the value of ref_t is 
relevant.

Mark


> 
> 
> 
> 
> 
> 
> 
> 
> 
> From: Justin A. Lemkul 
> To: Discussion list for GROMACS users 
> Sent: Tue, August 2, 2011 8:40:21 PM
> Subject: Re: [gmx-users] Problem with temperature and stochastic dynamics in 
> FEP
> 
> 
> 
> 
> Sanku M wrote:
> > Hi Justin,
> >   So, do you suggest that after minimization, I should generate the 
> >velocity at 300K instead ?
> > i.,e for equilibration,should following the set-up?
> > 
> 
> You should always equilibrate under the desired conditions.  Never make 
> sudden changes, or else you're going to damage the stability of the system.
> 
> -Justin
> 
> > gen-vel                  = yes
> >> >  gen-temp                = 300
> >> >  gen-seed                = -1
> > 
> > 
> > *From:* Justin A. Lemkul 
> > *To:* Discussion list for GROMACS users
>  
> > *Sent:* Tue, August 2, 2011 8:26:16 PM
> > *Subject:* Re: [gmx-users] Problem with temperature and stochastic dynamics 
> > in FEP
> > 
> > 
> > 
> > Justin A. Lemkul wrote:
> >  >
> >  >
> >  > Sanku M wrote:
> >  >> Hi Justin,
> >  >>  I first performed two minimization using steep and l-bfgs as the 
> >method on the solvated system for each lambda. Then I ran a 500ps NVT 
> >simulation on minimized system ( for each lambda) where essentially same 
> >.mdp file was used except the following change:
> >  >>    gen-vel                  = yes
> >  >>  gen-temp                = 280
> > 
>  >>  gen-seed     
>            = -1
> >  >>
> >  >> So, essentially, I heated the minimized system up from 280 K with ref_t 
> >300K using sd integrator. But, here also I found the average temperature 
> >goes to 303 K( instead of 300 K) .
> >  >
> >  > Unless you did simulated annealing in between, you didn't heat the 
> >system, you suddenly jolted the temperature, which is probably not stable.  
> >Maybe the integrator/Langevin thermostat does not react well to these 
> >conditions. Equilibrate under the conditions you wish to collect data.  Your 
> >data collection uses NPT, so running NVT at a different temperature and 
> >suddenly changing the temperature and introducing pressure coupling is not 
> >an appropriate procedure. If nothing else, your data will be unreliable, 
> >even if your temperature was correct.
> >  >
> > 
> > Sorry, misread the previous .mdp settings.  You're still using
>  NVT, but the point about the temperature is still valid here.
> > 
> > -Justin
> > 
> > -- 
> > 
> > Justin A. Lemkul
> > Ph.D. Candidate
> > ICTAS Doctoral Scholar
> > MILES-IGERT Trainee
> > Department of Biochemistry
> > Virginia Tech
> > Blacksburg, VA
> > jalemkul[at]vt.edu | (540) 231-9080
> > http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
> > 
> > 
> > -- gmx-users mailing list    gmx-users@gromacs.org 
> > 
> > http://lists.gromacs.org/mailman/listinfo/gmx-users
> > Please search the archive at 
> > http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> > Please don't post (un)subscribe requests to the list. Use the www interface 
> > or send it to gmx-users-requ...@gromacs.org 
> > .
> > Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> 
> --
>  
> 
> Justin A. Lemkul
> Ph.D. Candidate
> ICTAS Doctoral Scholar
> MILES-IGERT Trainee
> Department of Biochemistry
> Virginia
>  Tech
> Blacksburg, VA
> jalemkul[at]vt.edu | (540) 231-9080
> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
> 
> 
> -- gmx-users mailing list    gmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at 
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> Please don't post (un)subscribe requests to the list. Use the www inter

Re: [gmx-users] Problem with temperature and stochastic dynamics in FEP

[Sanku M]

Hi Justin,
   I did calculate the average temperature using last 2 ns of the 8ns data. It 
is still same 303K. 





From: Justin A. Lemkul 
To: Discussion list for GROMACS users 
Sent: Tue, August 2, 2011 9:55:40 PM
Subject: Re: [gmx-users] Problem with temperature and stochastic dynamics in FEP



Sanku M wrote:
> OKK, But, before, I have tried using md as integrator where I start the 
> initial 
>velocity as 280K ( which makes sense as the minimized configuration is 
>essentially at low temperature) and have used Nose-hoover at 300K and I never 
>had problem in getting the desired 300K after running a short NVT simulation 
>as 
>equilabration.
> The problem arises when sd integrator is used.
> In fact, I  just tried to generate the initial velocity at 300K along with sd 
>integrator. I got the 303 K as average temperature again. The same simulation 
>using md integrator always gives rise to 300K as average temperature( when 
>used 
>in conjunction with Nose-Hoover ) , even if I generate velocity at 280 K or 
>300K 
>.
> 
> I guess that the initial velocity-generation temperature should not be an 
> issue 
>here. The sd integrator in contrast to md integrator may be the problem.
> 
> 

I have never had such a problem when using sd.  I can always achieve the target 
temperature within 0.5 K, if not exactly, so I do not believe there is an 
inherent problem with sd.  Your first post indicated that your statistics were 
collected over the entire timeframe, which will not be correct for reasons I 
have already described.  The initial frames will represent unequilibrated data 
that may skew the averages.  What if you analyze only the last half of the 
time?  The last 1 ns?

-Justin

-- 

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin


-- gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search 
before posting!
Please don't post (un)subscribe requests to the list. Use the www interface or 
send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Problem with temperature and stochastic dynamics in FEP

[Justin A. Lemkul]

Sanku M wrote:

Hi Justin,
   I did calculate the average temperature using last 2 ns of the 8ns 
data. It is still same 303K. 





Going back to the original .mdp file, I can see that you're using some incorrect 
settings.  You've set rlist=rvdw even though you're using a shifted potential 
for van der Waals interactions.  Surely grompp should have warned you that this 
would lead to poor energy conservation.  Your solvent is probably heating up as 
a result and leading to the incorrect temperature.


-Justin



*From:* Justin A. Lemkul 
*To:* Discussion list for GROMACS users 
*Sent:* Tue, August 2, 2011 9:55:40 PM
*Subject:* Re: [gmx-users] Problem with temperature and stochastic 
dynamics in FEP




Sanku M wrote:
 > OKK, But, before, I have tried using md as integrator where I start 
the initial velocity as 280K ( which makes sense as the minimized 
configuration is essentially at low temperature) and have used 
Nose-hoover at 300K and I never had problem in getting the desired 300K 
after running a short NVT simulation as equilabration.

 > The problem arises when sd integrator is used.
 > In fact, I  just tried to generate the initial velocity at 300K along 
with sd integrator. I got the 303 K as average temperature again. The 
same simulation using md integrator always gives rise to 300K as average 
temperature( when used in conjunction with Nose-Hoover ) , even if I 
generate velocity at 280 K or 300K .

 >
 > I guess that the initial velocity-generation temperature should not 
be an issue here. The sd integrator in contrast to md integrator may be 
the problem.

 >
 >

I have never had such a problem when using sd.  I can always achieve the 
target temperature within 0.5 K, if not exactly, so I do not believe 
there is an inherent problem with sd.  Your first post indicated that 
your statistics were collected over the entire timeframe, which will not 
be correct for reasons I have already described.  The initial frames 
will represent unequilibrated data that may skew the averages.  What if 
you analyze only the last half of the time?  The last 1 ns?


-Justin

-- 

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin


-- gmx-users mailing listgmx-users@gromacs.org 


http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the www 
interface or send it to gmx-users-requ...@gromacs.org 
.

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


--


Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.

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

[gmx-users] Force-field checking options

[mcgrath]

Hi.  I'm fairly new to GROMACS, and I've been using it to run some
classical MD simulations.  In order to be sure that I'm using the
correct force field (I had to add a molecule to CHARMM27), I'm comparing
it to another simulation code that I know well, CP2K (using GROMACS
4.5.4 and a recent CVS version of CP2K).  Sadly, they are giving me
energy differences of about a factor of 2 for a 75000 atom protein+water
system.  As far as I can tell, I'm using the same PME parameters, and
there's not a big change in energy when I change those, anyway.  I've
been able to confirm that it's not a global problem (computing the
energy of only the 25,000 TIP3P waters give a result to within 1%, which
is not perfect, but better...I seem the same 1% difference if I only use
29 waters).  For a smaller system (using the molecule I added), the
total energy is incorrect, but the torsion and improper torsions are
good to within 1%.  So it looks like my topology is perhaps being
incorrectly specified, or the parameters for them.

What I would like to do is get GROMACS to print out all the charges (the
electrostatics/nonbonded are also different) that it is actually using,
as well as the force field parameters being used.  By using mdrun -v
-debug I get some of that, but lines like

c6= 1.48497790e-03, c12= 6.58807176e-06
c6= 3.78914556e-04, c12= 4.08982345e-07
c6= 2.02168059e-03, c12= 4.42785949e-06
c6= 1.71635114e-03, c12= 4.54480232e-06
c6= 2.60732602e-03, c12= 6.42257237e-06
c6= 3.75109637e-04, c12= 2.77185705e-07
c6= 1.98187144e-03, c12= 5.24788538e-06
c6= 6.18919948e-05, c12= 7.54611396e-09
c6= 1.73354731e-03, c12= 5.36550624e-06
c6= 4.41942073e-04, c12= 4.93156790e-07
c6= 6.79507386e-03, c12= 2.55060841e-05
c6= 1.61714898e-03, c12= 3.18964840e-06
c6= 2.48678902e-04, c12= 2.13336705e-07

are somewhat unclear.  They are obviously non-bonded terms, but
corresponding to which atoms?  The same goes for the bond angles and
torsions printed later.  The worst-case scenario is that I have to poke
around in the source files, which I would to avoid so I'm hoping there
is some documentation or more switches I can flip somewhere.  Thanks!

   Cheers, Matt

-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Force-field checking options

[Justin A. Lemkul]

mcgrath wrote:

Hi.  I'm fairly new to GROMACS, and I've been using it to run some
classical MD simulations.  In order to be sure that I'm using the
correct force field (I had to add a molecule to CHARMM27), I'm comparing
it to another simulation code that I know well, CP2K (using GROMACS
4.5.4 and a recent CVS version of CP2K).  Sadly, they are giving me
energy differences of about a factor of 2 for a 75000 atom protein+water
system.  As far as I can tell, I'm using the same PME parameters, and
there's not a big change in energy when I change those, anyway.  I've
been able to confirm that it's not a global problem (computing the
energy of only the 25,000 TIP3P waters give a result to within 1%, which
is not perfect, but better...I seem the same 1% difference if I only use
29 waters).  For a smaller system (using the molecule I added), the
total energy is incorrect, but the torsion and improper torsions are
good to within 1%.  So it looks like my topology is perhaps being
incorrectly specified, or the parameters for them.


Which energy terms show the discrepancy in the case of the twofold difference?



What I would like to do is get GROMACS to print out all the charges (the
electrostatics/nonbonded are also different) that it is actually using,
as well as the force field parameters being used.  By using mdrun -v
-debug I get some of that, but lines like



You don't need mdrun -debug for this.  The topology is static, so run gmxdump on 
your .tpr file.  This will map all parameters to each atom.


-Justin


c6= 1.48497790e-03, c12= 6.58807176e-06
c6= 3.78914556e-04, c12= 4.08982345e-07
c6= 2.02168059e-03, c12= 4.42785949e-06
c6= 1.71635114e-03, c12= 4.54480232e-06
c6= 2.60732602e-03, c12= 6.42257237e-06
c6= 3.75109637e-04, c12= 2.77185705e-07
c6= 1.98187144e-03, c12= 5.24788538e-06
c6= 6.18919948e-05, c12= 7.54611396e-09
c6= 1.73354731e-03, c12= 5.36550624e-06
c6= 4.41942073e-04, c12= 4.93156790e-07
c6= 6.79507386e-03, c12= 2.55060841e-05
c6= 1.61714898e-03, c12= 3.18964840e-06
c6= 2.48678902e-04, c12= 2.13336705e-07

are somewhat unclear.  They are obviously non-bonded terms, but
corresponding to which atoms?  The same goes for the bond angles and
torsions printed later.  The worst-case scenario is that I have to poke
around in the source files, which I would to avoid so I'm hoping there
is some documentation or more switches I can flip somewhere.  Thanks!

   Cheers, Matt



--


Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.

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

Re: [gmx-users] Problem with temperature and stochastic dynamics in FEP

[Sanku M]

Finally, I achieved the desired 300K average temperature using sd 
integrator 
. I had to reduce the time step from 0.002 ps to 0.001 ps and then running the 
NVT could produce 300K as average temperature.
Still Not sure why md integrator still could 300K average temperature using 
0.002 ps as time-step.




From: Justin A. Lemkul 
To: Discussion list for GROMACS users 
Sent: Tue, August 2, 2011 10:16:17 PM
Subject: Re: [gmx-users] Problem with temperature and stochastic dynamics in FEP



Sanku M wrote:
> Hi Justin,
>I did calculate the average temperature using last 2 ns of the 8ns data. 
> It 
>is still same 303K. 
>
> 

Going back to the original .mdp file, I can see that you're using some 
incorrect 
settings.  You've set rlist=rvdw even though you're using a shifted potential 
for van der Waals interactions.  Surely grompp should have warned you that this 
would lead to poor energy conservation.  Your solvent is probably heating up as 
a result and leading to the incorrect temperature.

-Justin

> 
> *From:* Justin A. Lemkul 
> *To:* Discussion list for GROMACS users 
> *Sent:* Tue, August 2, 2011 9:55:40 PM
> *Subject:* Re: [gmx-users] Problem with temperature and stochastic dynamics 
> in 
>FEP
> 
> 
> 
> Sanku M wrote:
>  > OKK, But, before, I have tried using md as integrator where I start the 
>initial velocity as 280K ( which makes sense as the minimized configuration is 
>essentially at low temperature) and have used Nose-hoover at 300K and I never 
>had problem in getting the desired 300K after running a short NVT simulation 
>as 
>equilabration.
>  > The problem arises when sd integrator is used.
>  > In fact, I  just tried to generate the initial velocity at 300K along with 
>sd integrator. I got the 303 K as average temperature again. The same 
>simulation 
>using md integrator always gives rise to 300K as average temperature( when 
>used 
>in conjunction with Nose-Hoover ) , even if I generate velocity at 280 K or 
>300K 
>.
>  >
>  > I guess that the initial velocity-generation temperature should not be an 
>issue here. The sd integrator in contrast to md integrator may be the problem.
>  >
>  >
> 
> I have never had such a problem when using sd.  I can always achieve the 
> target 
>temperature within 0.5 K, if not exactly, so I do not believe there is an 
>inherent problem with sd.  Your first post indicated that your statistics were 
>collected over the entire timeframe, which will not be correct for reasons I 
>have already described.  The initial frames will represent unequilibrated data 
>that may skew the averages.  What if you analyze only the last half of the 
>time?  The last 1 ns?
> 
> -Justin
> 
> -- 
> 
> Justin A. Lemkul
> Ph.D. Candidate
> ICTAS Doctoral Scholar
> MILES-IGERT Trainee
> Department of Biochemistry
> Virginia Tech
> Blacksburg, VA
> jalemkul[at]vt.edu | (540) 231-9080
> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
> 
> 
> -- gmx-users mailing listgmx-users@gromacs.org 
>
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at 
>http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> Please don't post (un)subscribe requests to the list. Use the www interface 
> or 
>send it to gmx-users-requ...@gromacs.org 
>.
> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

-- 

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin


-- gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search 
before posting!
Please don't post (un)subscribe requests to the list. Use the www interface or 
send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Blacksburg, VA
> jalemkul[at]vt.edu | (540) 231-9080
> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
> 
> 
> -- gmx-users mailing listgmx-users@gromacs.org 
>
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at 
>http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> Please don't post (un)subscribe requests to the list. Use the www interface 
> or 
>send it to gmx-users-requ...@gromacs.org 
>.
> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

-- 

Justin A. Lemkul
Ph.D.

Re: [gmx-users] Problem with temperature and stochastic dynamics in FEP

[Mark Abraham]

On 3/08/2011 2:23 PM, Sanku M wrote:


Finally, I achieved the desired 300K average temperature using sd 
integrator . I had to reduce the time step from 0.002 ps to 0.001 ps 
and then running the NVT could produce 300K as average temperature.
Still Not sure why md integrator still could 300K average temperature 
using 0.002 ps as time-step.


So your integration step was too large to conserve energy. There have 
been two suggestions made about why that might be true. Have you 
investigated them?


Mark




*From:* Justin A. Lemkul 
*To:* Discussion list for GROMACS users 
*Sent:* Tue, August 2, 2011 10:16:17 PM
*Subject:* Re: [gmx-users] Problem with temperature and stochastic 
dynamics in FEP




Sanku M wrote:
> Hi Justin,
>I did calculate the average temperature using last 2 ns of the 
8ns data. It is still same 303K.

>

Going back to the original .mdp file, I can see that you're using some 
incorrect settings.  You've set rlist=rvdw even though you're using a 
shifted potential for van der Waals interactions.  Surely grompp 
should have warned you that this would lead to poor energy 
conservation.  Your solvent is probably heating up as a result and 
leading to the incorrect temperature.


-Justin

> 
> *From:* Justin A. Lemkul mailto:jalem...@vt.edu>>
> *To:* Discussion list for GROMACS users >

> *Sent:* Tue, August 2, 2011 9:55:40 PM
> *Subject:* Re: [gmx-users] Problem with temperature and stochastic 
dynamics in FEP

>
>
>
> Sanku M wrote:
> > OKK, But, before, I have tried using md as integrator where I 
start the initial velocity as 280K ( which makes sense as the 
minimized configuration is essentially at low temperature) and have 
used Nose-hoover at 300K and I never had problem in getting the 
desired 300K after running a short NVT simulation as equilabration.

> > The problem arises when sd integrator is used.
> > In fact, I  just tried to generate the initial velocity at 300K 
along with sd integrator. I got the 303 K as average temperature 
again. The same simulation using md integrator always gives rise to 
300K as average temperature( when used in conjunction with Nose-Hoover 
) , even if I generate velocity at 280 K or 300K .

> >
> > I guess that the initial velocity-generation temperature should 
not be an issue here. The sd integrator in contrast to md integrator 
may be the problem.

> >
> >
>
> I have never had such a problem when using sd.  I can always achieve 
the target temperature within 0.5 K, if not exactly, so I do not 
believe there is an inherent problem with sd.  Your first post 
indicated that your statistics were collected over the entire 
timeframe, which will not be correct for reasons I have already 
described.  The initial frames will represent unequilibrated data that 
may skew the averages.  What if you analyze only the last half of the 
time?  The last 1 ns?

>
> -Justin
>
> -- 
>
> Justin A. Lemkul
> Ph.D. Candidate
> ICTAS Doctoral Scholar
> MILES-IGERT Trainee
> Department of Biochemistry
> Virginia Tech
> Blacksburg, VA
> jalemkul[at]vt.edu | (540) 231-9080
> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
>
> 
> -- gmx-users mailing list gmx-users@gromacs.org 
 >

> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> Please don't post (un)subscribe requests to the list. Use the www 
interface or send it to gmx-users-requ...@gromacs.org 
 
>.

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

-- 

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin


-- gmx-users mailing list gmx-users@gromacs.org 


http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the www 
interface or send it to gmx-users-requ...@gromacs.org 
.

Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Blacksburg, VA
> jalemkul[at]vt.edu | (540) 231-9080
> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
>
> 
> -- gmx-users mailing list gmx-users@gromacs.org 

[gmx-users] Problem of Kinetic Energy

[Size Zheng]

Hi, there,

I simulate a simple system of pure water with rigid TIP4P model in NVT at 273K.

Why the total kinetic energy, reported from g_energy, is not equal to
the sum of the rotational kinetic energy and the translational kinetic
energy, which are obtained from g_traj with options of -ekt and -ekr.

Even though I normalize the total kinetic energy (reported from g_energy)
with molecular number, its value is still not equal to the sum of the
rotational and translational kinetic energy ,obtained from g_traj
with options of -ekt and -ekr.

Many thanks,

Size Zheng--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Force-field checking options

[mcgrath]

Hi Justin, thanks for the response.  For the big system, the dihederal
is still fine, while now the improper, U-B, and nonbonded (CP2K groups
all nonbonded, 1-4, and real space Coulomb into the same term) are off
by a factor of 2-3.  Is there any way to print the energy due to angles
in GROMACS?  The nonbonded interactions dominate (for GROMACS, the short
range Coulomb term is an order of magnitude bigger than everything
else).  I guess it's possible my Ewald parameters are off, though that
wouldn't affect the SR term.  My box is 85x85x105 and I use a grid of
88x88x108, along with a 4th order interpolation for both simulations,
along with an ewald_rtol of 1e-5.  Playing with these doesn't seem to
change the numbers much, nor does switching between PME, SPME, or
regular Ewald (for CP2K).  rlist=rvdw=rcoulomb=1.2.  Is there a
parameter that I'm missing?  I've tried playing around with vdwtype and
coulombtype (I believe CP2K truncates and shifts nonbond and Coulombic
realspace potentials), but nothing had a significant effect.

Thanks for the gmxdump tip.  I'll start looking at the parameters for
the smaller system, and hopefully I can hunt down where the differences
are popping up. CP2K doesn't have CMAP implemented, but that only a
small part according to GROMACS.

  Cheers, Matt

> >  Forwarded Message 
> > From: Justin A. Lemkul 
> > Reply-to: jalem...@vt.edu, Discussion list for GROMACS users
> > 
> > To: Discussion list for GROMACS users 
> > Subject: Re: [gmx-users] Force-field checking options
> > Date: Wed, 03 Aug 2011 00:12:35 -0400
> > 
> > 
> > mcgrath wrote:
> > > Hi.  I'm fairly new to GROMACS, and I've been using it to run some
> > > classical MD simulations.  In order to be sure that I'm using the
> > > correct force field (I had to add a molecule to CHARMM27), I'm comparing
> > > it to another simulation code that I know well, CP2K (using GROMACS
> > > 4.5.4 and a recent CVS version of CP2K).  Sadly, they are giving me
> > > energy differences of about a factor of 2 for a 75000 atom protein+water
> > > system.  As far as I can tell, I'm using the same PME parameters, and
> > > there's not a big change in energy when I change those, anyway.  I've
> > > been able to confirm that it's not a global problem (computing the
> > > energy of only the 25,000 TIP3P waters give a result to within 1%, which
> > > is not perfect, but better...I seem the same 1% difference if I only use
> > > 29 waters).  For a smaller system (using the molecule I added), the
> > > total energy is incorrect, but the torsion and improper torsions are
> > > good to within 1%.  So it looks like my topology is perhaps being
> > > incorrectly specified, or the parameters for them.
> > 
> > Which energy terms show the discrepancy in the case of the twofold 
> > difference?
> > 
> > > 
> > > What I would like to do is get GROMACS to print out all the charges (the
> > > electrostatics/nonbonded are also different) that it is actually using,
> > > as well as the force field parameters being used.  By using mdrun -v
> > > -debug I get some of that, but lines like
> > > 
> > 
> > You don't need mdrun -debug for this.  The topology is static, so run 
> > gmxdump on 
> > your .tpr file.  This will map all parameters to each atom.
> > 
> > -Justin
> > 
> > > c6= 1.48497790e-03, c12= 6.58807176e-06
> > > c6= 3.78914556e-04, c12= 4.08982345e-07
> > > c6= 2.02168059e-03, c12= 4.42785949e-06
> > > c6= 1.71635114e-03, c12= 4.54480232e-06
> > > c6= 2.60732602e-03, c12= 6.42257237e-06
> > > c6= 3.75109637e-04, c12= 2.77185705e-07
> > > c6= 1.98187144e-03, c12= 5.24788538e-06
> > > c6= 6.18919948e-05, c12= 7.54611396e-09
> > > c6= 1.73354731e-03, c12= 5.36550624e-06
> > > c6= 4.41942073e-04, c12= 4.93156790e-07
> > > c6= 6.79507386e-03, c12= 2.55060841e-05
> > > c6= 1.61714898e-03, c12= 3.18964840e-06
> > > c6= 2.48678902e-04, c12= 2.13336705e-07
> > > 
> > > are somewhat unclear.  They are obviously non-bonded terms, but
> > > corresponding to which atoms?  The same goes for the bond angles and
> > > torsions printed later.  The worst-case scenario is that I have to poke
> > > around in the source files, which I would to avoid so I'm hoping there
> > > is some documentation or more switches I can flip somewhere.  Thanks!
> > > 
> > >Cheers, Matt
> > > 
> > 
> 

-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Force-field checking options

[Mark Abraham]

On 3/08/2011 3:34 PM, mcgrath wrote:

Hi Justin, thanks for the response.  For the big system, the dihederal
is still fine, while now the improper, U-B, and nonbonded (CP2K groups
all nonbonded, 1-4, and real space Coulomb into the same term) are off
by a factor of 2-3.  Is there any way to print the energy due to angles
in GROMACS?


Uh, they're available wherever you were getting energy break-downs from 
(e.g. in the .log file, or via g_energy on the .edr file). You can't get 
a break-down for each interaction, however. Simplify your system to 
probe things here. Do zero-step MD (not EM) without constraints, to 
evaluate the energy of a single conformation, and compare that with your 
other software. Complex things are complex to compare. :-) Reduce the 
complexity.



   The nonbonded interactions dominate (for GROMACS, the short
range Coulomb term is an order of magnitude bigger than everything
else).  I guess it's possible my Ewald parameters are off, though that
wouldn't affect the SR term.  My box is 85x85x105 and I use a grid of
88x88x108, along with a 4th order interpolation for both simulations,
along with an ewald_rtol of 1e-5.  Playing with these doesn't seem to
change the numbers much, nor does switching between PME, SPME, or
regular Ewald (for CP2K).  rlist=rvdw=rcoulomb=1.2.  Is there a
parameter that I'm missing?  I've tried playing around with vdwtype and
coulombtype (I believe CP2K truncates and shifts nonbond and Coulombic
realspace potentials), but nothing had a significant effect.

Thanks for the gmxdump tip.  I'll start looking at the parameters for
the smaller system, and hopefully I can hunt down where the differences
are popping up. CP2K doesn't have CMAP implemented, but that only a
small part according to GROMACS.


...and the magnitude of its energy contribution is available.

Mark



   Cheers, Matt


 Forwarded Message 
From: Justin A. Lemkul
Reply-to: jalem...@vt.edu, Discussion list for GROMACS users

To: Discussion list for GROMACS users
Subject: Re: [gmx-users] Force-field checking options
Date: Wed, 03 Aug 2011 00:12:35 -0400


mcgrath wrote:

Hi.  I'm fairly new to GROMACS, and I've been using it to run some
classical MD simulations.  In order to be sure that I'm using the
correct force field (I had to add a molecule to CHARMM27), I'm comparing
it to another simulation code that I know well, CP2K (using GROMACS
4.5.4 and a recent CVS version of CP2K).  Sadly, they are giving me
energy differences of about a factor of 2 for a 75000 atom protein+water
system.  As far as I can tell, I'm using the same PME parameters, and
there's not a big change in energy when I change those, anyway.  I've
been able to confirm that it's not a global problem (computing the
energy of only the 25,000 TIP3P waters give a result to within 1%, which
is not perfect, but better...I seem the same 1% difference if I only use
29 waters).  For a smaller system (using the molecule I added), the
total energy is incorrect, but the torsion and improper torsions are
good to within 1%.  So it looks like my topology is perhaps being
incorrectly specified, or the parameters for them.

Which energy terms show the discrepancy in the case of the twofold difference?


What I would like to do is get GROMACS to print out all the charges (the
electrostatics/nonbonded are also different) that it is actually using,
as well as the force field parameters being used.  By using mdrun -v
-debug I get some of that, but lines like


You don't need mdrun -debug for this.  The topology is static, so run gmxdump on
your .tpr file.  This will map all parameters to each atom.

-Justin


c6= 1.48497790e-03, c12= 6.58807176e-06
c6= 3.78914556e-04, c12= 4.08982345e-07
c6= 2.02168059e-03, c12= 4.42785949e-06
c6= 1.71635114e-03, c12= 4.54480232e-06
c6= 2.60732602e-03, c12= 6.42257237e-06
c6= 3.75109637e-04, c12= 2.77185705e-07
c6= 1.98187144e-03, c12= 5.24788538e-06
c6= 6.18919948e-05, c12= 7.54611396e-09
c6= 1.73354731e-03, c12= 5.36550624e-06
c6= 4.41942073e-04, c12= 4.93156790e-07
c6= 6.79507386e-03, c12= 2.55060841e-05
c6= 1.61714898e-03, c12= 3.18964840e-06
c6= 2.48678902e-04, c12= 2.13336705e-07

are somewhat unclear.  They are obviously non-bonded terms, but
corresponding to which atoms?  The same goes for the bond angles and
torsions printed later.  The worst-case scenario is that I have to poke
around in the source files, which I would to avoid so I'm hoping there
is some documentation or more switches I can flip somewhere.  Thanks!

Cheers, Matt



--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.

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

Re: [gmx-users] Problem of Kinetic Energy

[Mark Abraham]

On 3/08/2011 3:25 PM, Size Zheng wrote:

Hi, there,

I simulate a simple system of pure water with rigid TIP4P model in NVT at 273K.

Why the total kinetic energy, reported from g_energy, is not equal to
the sum of the rotational kinetic energy and the translational kinetic
energy, which are obtained from g_traj with options of -ekt and -ekr.


Read g_traj -h. The -ekt and -ekr options do not operate on atoms.



Even though I normalize the total kinetic energy (reported from g_energy)
with molecular number, its value is still not equal to the sum of the
rotational and translational kinetic energy ,obtained from g_traj
with options of -ekt and -ekr.


A stationary swarm of bees has high KE, and yet the swarm has low 
rotational and translational KE.


Mark
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.

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

Re: [gmx-users] Force-field checking options

[Mark Abraham]

On 3/08/2011 3:45 PM, Mark Abraham wrote:

On 3/08/2011 3:34 PM, mcgrath wrote:

Hi Justin, thanks for the response.  For the big system, the dihederal
is still fine, while now the improper, U-B, and nonbonded (CP2K groups
all nonbonded, 1-4, and real space Coulomb into the same term) are off
by a factor of 2-3.  Is there any way to print the energy due to angles
in GROMACS?


Uh, they're available wherever you were getting energy break-downs 
from (e.g. in the .log file, or via g_energy on the .edr file). You 
can't get a break-down for each interaction, however. Simplify your 
system to probe things here. Do zero-step MD (not EM) without 
constraints, to evaluate the energy of a single conformation, and 
compare that with your other software. Complex things are complex to 
compare. :-) Reduce the complexity.



   The nonbonded interactions dominate (for GROMACS, the short
range Coulomb term is an order of magnitude bigger than everything
else).  I guess it's possible my Ewald parameters are off, though that
wouldn't affect the SR term.  My box is 85x85x105 and I use a grid of
88x88x108, along with a 4th order interpolation for both simulations,
along with an ewald_rtol of 1e-5.  Playing with these doesn't seem to
change the numbers much, nor does switching between PME, SPME, or
regular Ewald (for CP2K).  rlist=rvdw=rcoulomb=1.2.  Is there a
parameter that I'm missing?  I've tried playing around with vdwtype and
coulombtype (I believe CP2K truncates and shifts nonbond and Coulombic
realspace potentials), but nothing had a significant effect.

Thanks for the gmxdump tip.  I'll start looking at the parameters for
the smaller system, and hopefully I can hunt down where the differences
are popping up. CP2K doesn't have CMAP implemented, but that only a
small part according to GROMACS.


...and the magnitude of its energy contribution is available.


...and trivial to turn off in pdb2gmx.

Mark


Mark



   Cheers, Matt


 Forwarded Message 
From: Justin A. Lemkul
Reply-to: jalem...@vt.edu, Discussion list for GROMACS users

To: Discussion list for GROMACS users
Subject: Re: [gmx-users] Force-field checking options
Date: Wed, 03 Aug 2011 00:12:35 -0400


mcgrath wrote:

Hi.  I'm fairly new to GROMACS, and I've been using it to run some
classical MD simulations.  In order to be sure that I'm using the
correct force field (I had to add a molecule to CHARMM27), I'm 
comparing

it to another simulation code that I know well, CP2K (using GROMACS
4.5.4 and a recent CVS version of CP2K).  Sadly, they are giving me
energy differences of about a factor of 2 for a 75000 atom 
protein+water

system.  As far as I can tell, I'm using the same PME parameters, and
there's not a big change in energy when I change those, anyway.  I've
been able to confirm that it's not a global problem (computing the
energy of only the 25,000 TIP3P waters give a result to within 1%, 
which
is not perfect, but better...I seem the same 1% difference if I 
only use

29 waters).  For a smaller system (using the molecule I added), the
total energy is incorrect, but the torsion and improper torsions are
good to within 1%.  So it looks like my topology is perhaps being
incorrectly specified, or the parameters for them.
Which energy terms show the discrepancy in the case of the twofold 
difference?


What I would like to do is get GROMACS to print out all the 
charges (the
electrostatics/nonbonded are also different) that it is actually 
using,

as well as the force field parameters being used.  By using mdrun -v
-debug I get some of that, but lines like

You don't need mdrun -debug for this.  The topology is static, so 
run gmxdump on

your .tpr file.  This will map all parameters to each atom.

-Justin


c6= 1.48497790e-03, c12= 6.58807176e-06
c6= 3.78914556e-04, c12= 4.08982345e-07
c6= 2.02168059e-03, c12= 4.42785949e-06
c6= 1.71635114e-03, c12= 4.54480232e-06
c6= 2.60732602e-03, c12= 6.42257237e-06
c6= 3.75109637e-04, c12= 2.77185705e-07
c6= 1.98187144e-03, c12= 5.24788538e-06
c6= 6.18919948e-05, c12= 7.54611396e-09
c6= 1.73354731e-03, c12= 5.36550624e-06
c6= 4.41942073e-04, c12= 4.93156790e-07
c6= 6.79507386e-03, c12= 2.55060841e-05
c6= 1.61714898e-03, c12= 3.18964840e-06
c6= 2.48678902e-04, c12= 2.13336705e-07

are somewhat unclear.  They are obviously non-bonded terms, but
corresponding to which atoms?  The same goes for the bond angles and
torsions printed later.  The worst-case scenario is that I have to 
poke
around in the source files, which I would to avoid so I'm hoping 
there

is some documentation or more switches I can flip somewhere.  Thanks!

Cheers, Matt





--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.