[gmx-users] Exploding temp/pressure.

[Nash, Anthony]

Hi all,

I¹m trying to equilibrate a Martini CG simulation from an initial
atomistic structure. Eq and Fc values were derived using an atomistic
system. I¹ve started the dt at 0.0005 for 60 steps, moving through
0.001, 0.0015 and 0.002 for the same number of steps, using the .mdp
details below (with the exception of the time step and the output
delimiters - the remaining settings are from the Martini website). During
these steps the temperature dropped from approx 6K down to 310K and
the potential energy dropped from positive figures down to around
-1.4*10^6. A final 60 steps with dt = 0.002 showed a converged temp,
pressure, volume, potential energy (and total energy), and very little
change in the unit dimensions.

I now want to move the time step up to 0.025. I¹ve tried various dt values
(0.005, 0.01, 0.015Š 0.025) using the below setup, but each time dt causes
almost immediately structural explosion. The temperature and the pressure
rockets up. I changed the pressure coupling to berendsen and the same
happens. I have been through dozens of tau_p and tau_t combinations, each
one yields a trajectory of varying length before exploding (the protein
Œexplodes¹ at the exact time the the pressure and temperature increases).
I have also changed the thermostat to berendsen and Nose-Hoover, using a
range of 5*(dt*nsttcouple) and 20(dt*nsttcouple), respectively.


I would appreciate any thoughts or suggestions on how to solve my
temperature/pressure problem when I try increasing the dt to a CG value.

Thanks
Anthony


title= Martini
integrator   = md
dt   = 0.025
nsteps   = 30
nstcomm  = 100


nstxout  = 1
nstvout  = 1
nstfout  = 0
nstlog   = 1
nstenergy= 1
nstxout-compressed   = 0
compressed-x-precision   = 0
;compressed-x-grps=
energygrps   = collagen solvent
cutoff-scheme= Verlet
nstlist  = 20
ns_type  = grid
pbc  = xyz
verlet-buffer-tolerance  = 0.005
coulombtype  = PME ;reaction-field
rcoulomb = 1.1
fourierspacing   = 0.16 ;0.2  ;0.12
epsilon_r= 15   ; 2.5 (with polarizable water)
epsilon_rf   = 0
vdw_type = cutoff
vdw-modifier = Potential-shift-verlet
rvdw = 1.1

tcoupl   = v-rescale ;berendsen ;v-rescale
tc-grps  = collagen solvent
tau_t= 0.4 0.4 ;2.0 2.0
ref_t= 310 310

Pcoupl   = parrinello-rahman
Pcoupltype   = semiisotropic
tau_p= 12.0  ;parrinello-rahman is more stable with
larger tau-p, DdJ, 20130422
compressibility  = 4.5e-5 0
ref_p= 1.0 1.0
refcoord_scaling = com

gen_vel  = no
gen_temp = 310
gen_seed = 473529
continuation = yes
constraints  = none
constraint_algorithm = lincs
lincs-warnangle = 45
lincs-order=18
lincs-iter=4



-- 
Gromacs Users mailing list

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

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

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

Re: [gmx-users] CG Lincs errors

[Nash, Anthony]

Hi Peter, 

Thanks for the reply, I know we spoke in length on this mater only just
recently. Many thanks for that.

I’ve taken the time step of collagen in vacu down to 0.0001 and I’ve
dropped the temp down to 280. I hope, running over 16 cores for two days
that this should relieve any tension in the backdone before in gradually
increase to 20-40fs. Again, all in vacu with no solvent.

I’ve actually thought about writing a script to modify the equilibrium
bond angles in the CG.itp file for the backbone, using the atomistic
structure as a template. After all, true collagen does not replicate
‘ideal’ collagen along the stretch of the protein e.g., the MMP1 binding
site is not very tightly coiled. Perhaps, there lies my problem.

I haven’t thought too much of the actual full fibril structure, I want to
capture the D-band gap in  type I collagen. When I give it some more
thought I will probably look into semi-isotropic pressure coupling first.

Thanks
Anthony 

Dr Anthony Nash
Department of Chemistry
University College London

Skeletal Tissue Dynamics Group
Committee member of London Matrix Group @LondonMatrixGrp





On 16/12/2016 09:02, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Peter Kroon" <gromacs.org_gmx-users-boun...@maillist.sys.kth.se
on behalf of p.c.kr...@rug.nl> wrote:

>As a note to Alex (and the rest of the list), the coarse-grained Martini
>forcefield is usually run with timesteps between 20-40 fs. 15fs is
>already rather low. I do agree that longer equilibration at low timestep
>(5 or 10) might help.
>
>Alternatively, Do you think a semiisotropic pressure coupling might be
>applicable in this case, since it's an infinite collagen polymer?
>
>
>Peter (PhD in the Martini group)
>
>
>On 16-12-16 00:21, Nash, Anthony wrote:
>> Alex and Mark, thanks for the information. I’ll drop dt down,
>> significantly, drop the temperature, and run it for a long time.
>>
>> Thanks for the ideas.
>> Anthony 
>>
>>
>> On 15/12/2016 21:54, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se
>>on
>> behalf of Alex" <gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
>> behalf of nedoma...@gmail.com> wrote:
>>
>>> Mark is right, no two ways about it. For initial equilibration and
>>> assessing preexisting structural strains try vacuum, _much_ smaller
>>> timesteps and possibly low temperatures in vacuum, only then transfer
>>>to
>>> solvent, etc. Algorithmically, LINCS requires convergence and you
>>>already
>>> are using a pretty high LINCS order... From what I see, dt = 15 fs at
>>>310K
>>> looks like a cowboy mode simulation in this case.
>>>
>>> Alex
>>>
>>> On Thu, Dec 15, 2016 at 2:32 PM, Mark Abraham
>>><mark.j.abra...@gmail.com>
>>> wrote:
>>>
>>>> Hi,
>>>>
>>>> If a simulation isn't stable with a small time step (as I think you
>>>>are
>>>> saying) then moving to a larger time step is guaranteed to make that
>>>> worse.
>>>> Try an even smaller time step, for a long time, and see what happens.
>>>>Or
>>>> take a subset of your protein and see what happens. Or simulate in
>>>>vacuo
>>>> for a while. Your topology could be unsuited to your starting
>>>>structure,
>>>> e.g. some part is under a lot of tension that gets released at some
>>>> point
>>>> and no finite time step can in practice deal with the velocity of the
>>>> recoil...
>>>>
>>>> Mark
>>>>
>>>> On Thu, 15 Dec 2016 23:06 Nash, Anthony <a.n...@ucl.ac.uk> wrote:
>>>>
>>>>> Hi all,
>>>>>
>>>>> I¹m hoping for some help. I¹m very sorry, this is a bit of a long
>>>>>one.
>>>>>
>>>>> I¹ve been struggling for almost a month trying to run a CG
>>>> representation
>>>>> of our all-atom model of a collagen protein (3 polypeptide chains in
>>>>>a
>>>>> protein). Our original AMBER all-atom model had been successful
>>>> modelling
>>>>> using MD. I went on to use the latest version of Martinize.py with
>>>>>the
>>>>> latest version of the MARTINI forcefield fields.
>>>>>
>>>>> After a little tweaking (the way AMBER names histidine residues), I
>>>>> successful converted the molecule (approx 3100 amino acids) into a CG
>>>>> representation. I successfully energy min the protein in vacuum to a
>>>>> threshold of 500, and i

Re: [gmx-users] CG Lincs errors

[Nash, Anthony]

Alex and Mark, thanks for the information. I’ll drop dt down,
significantly, drop the temperature, and run it for a long time.

Thanks for the ideas.
Anthony 


On 15/12/2016 21:54, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Alex" <gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of nedoma...@gmail.com> wrote:

>Mark is right, no two ways about it. For initial equilibration and
>assessing preexisting structural strains try vacuum, _much_ smaller
>timesteps and possibly low temperatures in vacuum, only then transfer to
>solvent, etc. Algorithmically, LINCS requires convergence and you already
>are using a pretty high LINCS order... From what I see, dt = 15 fs at 310K
>looks like a cowboy mode simulation in this case.
>
>Alex
>
>On Thu, Dec 15, 2016 at 2:32 PM, Mark Abraham <mark.j.abra...@gmail.com>
>wrote:
>
>> Hi,
>>
>> If a simulation isn't stable with a small time step (as I think you are
>> saying) then moving to a larger time step is guaranteed to make that
>>worse.
>> Try an even smaller time step, for a long time, and see what happens. Or
>> take a subset of your protein and see what happens. Or simulate in vacuo
>> for a while. Your topology could be unsuited to your starting structure,
>> e.g. some part is under a lot of tension that gets released at some
>>point
>> and no finite time step can in practice deal with the velocity of the
>> recoil...
>>
>> Mark
>>
>> On Thu, 15 Dec 2016 23:06 Nash, Anthony <a.n...@ucl.ac.uk> wrote:
>>
>> > Hi all,
>> >
>> > I¹m hoping for some help. I¹m very sorry, this is a bit of a long one.
>> >
>> > I¹ve been struggling for almost a month trying to run a CG
>>representation
>> > of our all-atom model of a collagen protein (3 polypeptide chains in a
>> > protein). Our original AMBER all-atom model had been successful
>>modelling
>> > using MD. I went on to use the latest version of Martinize.py with the
>> > latest version of the MARTINI forcefield fields.
>> >
>> > After a little tweaking (the way AMBER names histidine residues), I
>> > successful converted the molecule (approx 3100 amino acids) into a CG
>> > representation. I successfully energy min the protein in vacuum to a
>> > threshold of 500, and in solvent to a threshold of 750 using steepest
>> > descent. Looking for a system at an energy min of a threshold around
>>300
>> I
>> > begin to see LINCS warnings. Observing the initial structure, there is
>> > nothing obviously wrong with the bond network (both protein and
>>polarized
>> > CG water).
>> >
>> > I take the system that energy mins at 750 (protein-water mix, with no
>> > fault reported), and went straight to NPT, 20fs step. Blew up. After a
>> bit
>> > of chatting with the MARTINI community, I¹ve started with an NVT
>> ensemble,
>> > beginning at 5s then through 10fs, 15fs, and 20fs. I only run for 1000
>> > steps before switching. Keeping any of the simulations running for
>>longer
>> > throws lincs warnings followed by a segmentation fault from the
>>warning:
>> >
>> > "3 particles communicated to PME rank 7 are more than 2/3 times the
>> > cut-off out of the domain decomposition cell of their charge group in
>> > dimension x."
>> >
>> > Observing the trajectories of any of the extended simulations shows
>>the
>> > protein snapping like a rope, and always at the same place. I have
>> watched
>> > every trajectory at this point, using numerous energy min start
>>points,
>> to
>> > try and understand why it is blowing up. I can¹t see any obvious
>>reason.
>> I
>> > was told to consider how the temperature is changing. Below is an
>>example
>> > of the temperature and pressure from an NPT of 20fs step continued
>>from
>> > the very short 20fs step NVT simulation (hoping that perhaps CG
>>without
>> > pressure just doesn¹t behave happily; I was wrong).
>> >
>> >
>> > TEMP:
>> > Š
>> > 6.63  311.000336
>> > 6.645000  311.371643
>> > 6.66  311.724213
>> > 6.675000  313.878693
>> > 6.69  681558.937500
>> >
>> >
>> > PRESSURE:
>> > Š
>> > 6.633.559879
>> > 6.6450003.901433
>> > 6.663.589078
>> > 6.6750004.158611
>> > 6.69  81762.437500
>> >
>> > 

[gmx-users] CG Lincs errors

[Nash, Anthony]

Hi all,

I¹m hoping for some help. I¹m very sorry, this is a bit of a long one.

I¹ve been struggling for almost a month trying to run a CG representation
of our all-atom model of a collagen protein (3 polypeptide chains in a
protein). Our original AMBER all-atom model had been successful modelling
using MD. I went on to use the latest version of Martinize.py with the
latest version of the MARTINI forcefield fields.

After a little tweaking (the way AMBER names histidine residues), I
successful converted the molecule (approx 3100 amino acids) into a CG
representation. I successfully energy min the protein in vacuum to a
threshold of 500, and in solvent to a threshold of 750 using steepest
descent. Looking for a system at an energy min of a threshold around 300 I
begin to see LINCS warnings. Observing the initial structure, there is
nothing obviously wrong with the bond network (both protein and polarized
CG water).

I take the system that energy mins at 750 (protein-water mix, with no
fault reported), and went straight to NPT, 20fs step. Blew up. After a bit
of chatting with the MARTINI community, I¹ve started with an NVT ensemble,
beginning at 5s then through 10fs, 15fs, and 20fs. I only run for 1000
steps before switching. Keeping any of the simulations running for longer
throws lincs warnings followed by a segmentation fault from the warning:

"3 particles communicated to PME rank 7 are more than 2/3 times the
cut-off out of the domain decomposition cell of their charge group in
dimension x."

Observing the trajectories of any of the extended simulations shows the
protein snapping like a rope, and always at the same place. I have watched
every trajectory at this point, using numerous energy min start points, to
try and understand why it is blowing up. I can¹t see any obvious reason. I
was told to consider how the temperature is changing. Below is an example
of the temperature and pressure from an NPT of 20fs step continued from
the very short 20fs step NVT simulation (hoping that perhaps CG without
pressure just doesn¹t behave happily; I was wrong).


TEMP:
Š
6.63  311.000336
6.645000  311.371643
6.66  311.724213
6.675000  313.878693
6.69  681558.937500


PRESSURE:
Š
6.633.559879
6.6450003.901433
6.663.589078
6.6750004.158611
6.69  81762.437500

The final LINCS warning from this same run:

Step 300, time 4.5 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 0.35, max 0.003386 (between atoms 2125 and 2126)
bonds that rotated more than 45 degrees:
 atom 1 atom 2  angle  previous, current, constraint length
   2125   2126   68.30.2781   0.2691  0.2700
   2125   2127   45.90.2789   0.2701  0.2700


At this stage the structure ruptures as described above.


My NVT settings (with NPT included to save space) are:

-
title= Martini

integrator   = md
dt   = 0.015
nsteps   = 1000
nstcomm  = 100
;comm-grps   =

nstxout  = 1000
nstvout  = 1000
nstfout  = 0
nstlog   = 1
nstenergy= 1
nstxout-compressed   = 0
compressed-x-precision   = 0
;compressed-x-grps=
energygrps   = collagen solvent

cutoff-scheme= Verlet
nstlist  = 20
ns_type  = grid

pbc  = xyz
verlet-buffer-tolerance  = 0.005

coulombtype  = PME ;reaction-field
rcoulomb = 1.1
fourierspacing   = 0.16 ;0.2  ;0.12

epsilon_r= 2.5 ;15  ; 2.5 (with polarizable water)
epsilon_rf   = 0
vdw_type = cutoff
vdw-modifier = Potential-shift-verlet
rvdw = 1.1

tcoupl   = v-rescale ;berendsen ;v-rescale
tc-grps  = collagen solvent
tau_t= 0.5 0.5 ;1.0 1.0
ref_t= 310 310

Pcoupl   = berendsen   ;parrinello-rahman
Pcoupltype   = isotropic
tau_p= 12.0  ; parrinello-rahman is more stable with
larger tau-p, DdJ, 20130422
compressibility  = 10e-4
ref_p= 1.0
refcoord_scaling = com

gen_vel  = no
gen_temp = 310
gen_seed = 473529

continuation = yes
constraints  = none
constraint_algorithm = lincs
lincs-warnangle = 45
lincs-order=8
lincs-iter=4


‹‹

Every setting bar the lincs iter, order, warnangle were supplied with the
latest version of MARTINI. During many NVT runs I have adjusted the tau-t
to try and keep the thermostat from oscillating its way into infinity.

I¹m curious, will an out of control thermostat break a structure, or will
a structure breaking (for what ever reason this structure is breaking)
cause the thermostat to go out of control?

My only 

Re: [gmx-users] Fine tune the RDF of water around a dummy metal

[Nash, Anthony]

I've only been able to make contact with a colleague of the original author, 
and they have been extremely helpful with suggestions but the values still lies 
out range of that published. 

Therefore, I was hoping for some useful suggestion from individuals with more 
experience than I in the guts of Gromacs on what parameter beyond, charge, 
sigma, epsilon, mass, FC and eq distances, could contribute to a shift in RDF 
from metal to oxygen of water. 

Thanks
Anthony 

Dr Anthony Nash
Department of Chemistry
University College London

From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se 
[gromacs.org_gmx-users-boun...@maillist.sys.kth.se] on behalf of Justin Lemkul 
[jalem...@vt.edu]
Sent: 12 October 2016 22:56
To: gmx-us...@gromacs.org
Subject: Re: [gmx-users] Fine tune the RDF of water around a dummy metal

On 10/12/16 5:32 PM, Nash, Anthony wrote:
> Hi all,
>
> I¹m trying to fine tune the rdf of tip3p water molecules around a central
> metal dummy molecule ("Force Field Independent Metal Parameters Using a
> Nonbonded Dummy Model²), essentially a central metal (with vdw parameter
> and -1 charge) covalently bonded to six Œdummy¹ atoms (no vdw parameters,
> each with a +0.5 charge). Despite all my efforts at recreating the RDF
> from the published work, my first peak deviates by almost 0.8 Angstroms.
> The original publication showed that a deviation from their rdf resulted
> in a large deviation from the correct solvation free energy. Also, a metal
> demonstrating a deviation of around 0.8 angstrom in a water box will more
> than likely prevent the metal from settling in plane to the nitrogen atoms
> of a heme group when the parameters are put into practice.
>
> I have looked into adjusting the following in my attempt at bringing the
> first hydration shell closer to the dummy atom:
>
> -adjusting the sigma & epsilon (started from the published values I¹ve
> iteratively gone through tens of paired values, with sigma never getting
> below 0.48129. Any attempt at doing so causes steepest descent to crash
> during an energy minimisation)
> -shortened the equilibrium bond lengths of the central metal atom to its
> bonded dummy atoms
> -increased the positive charge on the dummy atoms whilst reducing the
> negative charge on the central metal atom
>
> 0.32 Angstrom is as close as I can get it without the energy minimisation
> from crashing out (excessive force).
>
> A bit of a long shot but any suggestions on how I could tune these
> parameters further?
>

What do the authors say when you tell them you can't reproduce their published
results?

-Justin

--
==

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

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

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

==
--
Gromacs Users mailing list

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

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

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

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

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

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

[gmx-users] Fine tune the RDF of water around a dummy metal

[Nash, Anthony]

Hi all,

I¹m trying to fine tune the rdf of tip3p water molecules around a central
metal dummy molecule ("Force Field Independent Metal Parameters Using a
Nonbonded Dummy Model²), essentially a central metal (with vdw parameter
and -1 charge) covalently bonded to six Œdummy¹ atoms (no vdw parameters,
each with a +0.5 charge). Despite all my efforts at recreating the RDF
from the published work, my first peak deviates by almost 0.8 Angstroms.
The original publication showed that a deviation from their rdf resulted
in a large deviation from the correct solvation free energy. Also, a metal
demonstrating a deviation of around 0.8 angstrom in a water box will more
than likely prevent the metal from settling in plane to the nitrogen atoms
of a heme group when the parameters are put into practice.

I have looked into adjusting the following in my attempt at bringing the
first hydration shell closer to the dummy atom:

-adjusting the sigma & epsilon (started from the published values I¹ve
iteratively gone through tens of paired values, with sigma never getting
below 0.48129. Any attempt at doing so causes steepest descent to crash
during an energy minimisation)
-shortened the equilibrium bond lengths of the central metal atom to its
bonded dummy atoms
-increased the positive charge on the dummy atoms whilst reducing the
negative charge on the central metal atom

0.32 Angstrom is as close as I can get it without the energy minimisation
from crashing out (excessive force).

A bit of a long shot but any suggestions on how I could tune these
parameters further?


Dr Anthony Nash
Department of Chemistry
University College London

Skeletal Tissue Dynamics Group
Committee member of London Matrix Group @LondonMatrixGrp

-- 
Gromacs Users mailing list

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

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

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

Re: [gmx-users] energy minimisation - LINCS WARNING

[Nash, Anthony]

That’s a good suggestion, Mark. Will do :-)

Thanks
Anthony

Dr Anthony Nash
Department of Chemistry
University College London

Skeletal Tissue Dynamics Group
Committee member of London Matrix Group @LondonMatrixGrp





On 03/10/2016 12:03, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Mark Abraham" <gromacs.org_gmx-users-boun...@maillist.sys.kth.se
on behalf of mark.j.abra...@gmail.com> wrote:

>Hi,
>
>Sounds like you could check the output from that tool, and if it gave no
>warning (or offers you no way to choose another rotamer) then you could
>make some constructive feedback to its authors :-)
>
>Mark
>
>On Mon, 3 Oct 2016 11:22 Nash, Anthony <a.n...@ucl.ac.uk> wrote:
>
>>
>> Hi Justin,
>>
>>
>> I’ve tried out all of your suggestions, they worked to an extent but the
>> in vacu idea was a good call as it’s actually helped me isolate what the
>> real problem is: I-TASSER has put the C-N backbone bond of two residues
>> through the ring of an adjacent PHE residue!
>>
>> Thanks
>> Anthony
>>
>>
>> On 02/10/2016 23:16, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se
>>on
>> behalf of Justin Lemkul"
>> <gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
>> jalem...@vt.edu> wrote:
>>
>> >
>> >Please don't reply to digests or unrelated threads when starting a new
>> >topic; it
>> >creates a mess in the archive.
>> >
>> >On 10/2/16 6:11 PM, Nash, Anthony wrote:
>> >> Hi all,
>> >>
>> >> I had a homology/de-novo model .pdb converted into .gro, solvated,
>> >> neutralised and now I¹m going through a series of energy minimisation
>> >> steps. Unfortunately, during energy minimisation I get LINCS WARNINGS
>> >> (angle relative constraint deviation). The the naked eye, the atoms
>> >> involved don¹t look untoward.
>> >>
>> >>
>> >> When I used LINCS all-bond I get the following error:
>> >> ---
>> >>
>> >> Step -1, time -0.001 (ps)  LINCS WARNING
>> >> relative constraint deviation after LINCS:
>> >> rms 0.001657, max 0.033500 (between atoms 1362 and 1363)
>> >> bonds that rotated more than 30 degrees:
>> >>  atom 1 atom 2  angle  previous, current, constraint length
>> >>
>> >>
>> >>
>> >> When I reduced this to only be the hydrogen bonds (h-bonds) I get:
>> >> ---
>> >> Step 20, time 0.02 (ps)  LINCS WARNING
>> >> relative constraint deviation after LINCS:
>> >> rms 0.00, max 0.00 (between atoms 1322 and 1324)
>> >> bonds that rotated more than 30 degrees:
>> >>  atom 1 atom 2  angle  previous, current, constraint length
>> >>1360   1361   37.50.1010   0.1010  0.1010
>> >> [continues with more errors until Segmentation fault 11]
>> >>
>> >>
>> >> I have tried both constraints with a reduced and increased deltaStep
>> >> (0.2,0.02,0.002,0.0002) no change. I have also taken constraints off,
>> >>this
>> >> does work but then fails again when constraints are returned. I have
>> >>tried
>> >> both Steep and cg, both fail with the same result.
>> >>
>> >> The .mdp file I used (taken from the latest release of Justin¹s
>> >>tutorials)
>> >> is:
>> >> ; energy.mdp - used as input into grompp to generate em.tpr
>> >> ; Parameters describing what to do, when to stop and what to save
>> >> integrator  = steep ; Algorithm (steep = steepest descent
>> >> minimization)
>> >> emtol   = 700.0; Stop minimization when the maximum
>> >>force
>> >> < 1000.0 kJ/mol/nm
>> >> emstep  = 0.002  ; Energy step size
>> >> nsteps  = 20; Maximum number of
>> >>(minimization)
>> >> steps to perform
>> >>
>> >> ; Parameters describing how to find the neighbors of each atom and
>>how
>> >>to
>> >> calculate the interactions
>> >> cutoff_scheme = verlet
>> >> nstlist = 20; Frequency to update the neighbor
>>list
>> >> and long range forces
>> >> ns_type = grid  ; Method to determine neighbor list
>> >> (simple, grid)
>> >> rlist   = 1

Re: [gmx-users] energy minimisation - LINCS WARNING

[Nash, Anthony]

Hi Justin,


I’ve tried out all of your suggestions, they worked to an extent but the
in vacu idea was a good call as it’s actually helped me isolate what the
real problem is: I-TASSER has put the C-N backbone bond of two residues
through the ring of an adjacent PHE residue!

Thanks
Anthony 


On 02/10/2016 23:16, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Justin Lemkul"
<gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
jalem...@vt.edu> wrote:

>
>Please don't reply to digests or unrelated threads when starting a new
>topic; it 
>creates a mess in the archive.
>
>On 10/2/16 6:11 PM, Nash, Anthony wrote:
>> Hi all,
>>
>> I had a homology/de-novo model .pdb converted into .gro, solvated,
>> neutralised and now I¹m going through a series of energy minimisation
>> steps. Unfortunately, during energy minimisation I get LINCS WARNINGS
>> (angle relative constraint deviation). The the naked eye, the atoms
>> involved don¹t look untoward.
>>
>>
>> When I used LINCS all-bond I get the following error:
>> ---
>>
>> Step -1, time -0.001 (ps)  LINCS WARNING
>> relative constraint deviation after LINCS:
>> rms 0.001657, max 0.033500 (between atoms 1362 and 1363)
>> bonds that rotated more than 30 degrees:
>>  atom 1 atom 2  angle  previous, current, constraint length
>>
>>
>>
>> When I reduced this to only be the hydrogen bonds (h-bonds) I get:
>> ---
>> Step 20, time 0.02 (ps)  LINCS WARNING
>> relative constraint deviation after LINCS:
>> rms 0.00, max 0.00 (between atoms 1322 and 1324)
>> bonds that rotated more than 30 degrees:
>>  atom 1 atom 2  angle  previous, current, constraint length
>>1360   1361   37.50.1010   0.1010  0.1010
>> [continues with more errors until Segmentation fault 11]
>>
>>
>> I have tried both constraints with a reduced and increased deltaStep
>> (0.2,0.02,0.002,0.0002) no change. I have also taken constraints off,
>>this
>> does work but then fails again when constraints are returned. I have
>>tried
>> both Steep and cg, both fail with the same result.
>>
>> The .mdp file I used (taken from the latest release of Justin¹s
>>tutorials)
>> is:
>> ; energy.mdp - used as input into grompp to generate em.tpr
>> ; Parameters describing what to do, when to stop and what to save
>> integrator  = steep ; Algorithm (steep = steepest descent
>> minimization)
>> emtol   = 700.0; Stop minimization when the maximum
>>force
>> < 1000.0 kJ/mol/nm
>> emstep  = 0.002  ; Energy step size
>> nsteps  = 20; Maximum number of
>>(minimization)
>> steps to perform
>>
>> ; Parameters describing how to find the neighbors of each atom and how
>>to
>> calculate the interactions
>> cutoff_scheme = verlet
>> nstlist = 20; Frequency to update the neighbor list
>> and long range forces
>> ns_type = grid  ; Method to determine neighbor list
>> (simple, grid)
>> rlist   = 1.4   ; Cut-off for making neighbor list
>>(short
>> range forces)
>> coulombtype = PME   ; Treatment of long range electrostatic
>> interactions
>> rcoulomb= 1.4   ; Short-range electrostatic cut-off
>> rvdw= 1.4   ; Short-range Van der Waals cut-off
>> pbc = xyz   ; Periodic Boundary Conditions (yes/no)
>>
>> ;ADDED THIS BIT MYSELF
>> continuation= no   ; first dynamics run
>> constraint_algorithm = lincs; holonomic constraints
>> constraints = h-bonds; all-bonds (even heavy atom-H bonds)
>>constrained
>> lincs_iter  = 1 ; accuracy of LINCS
>> lincs_order = 8 ; also related to accuracy
>>
>>
>> The .pdb was generated using I-TASSER. I have tried refining the
>>hydrogen
>> placement via -ignh during pdb2gmx, and I also tried ³Protein
>>Preparation²
>> (bond ordering - again, hydrogen bond placement) in Schrodinger. All the
>> same.
>>
>> I was hoping that I could actually watch a rough-and-ready minimisation
>> using Avogadro, however, the release I have keeps crashing whilst trying
>> to open large structures.
>>
>> Any suggestions on what I could try next?
>>
>
>All signs point to unreasonable starting structure.  How about minimizing
>the 
>protein alone, in vacuo?  Maybe re

[gmx-users] energy minimisation - LINCS WARNING

[Nash, Anthony]

Hi all,

I had a homology/de-novo model .pdb converted into .gro, solvated,
neutralised and now I¹m going through a series of energy minimisation
steps. Unfortunately, during energy minimisation I get LINCS WARNINGS
(angle relative constraint deviation). The the naked eye, the atoms
involved don¹t look untoward.


When I used LINCS all-bond I get the following error:
---

Step -1, time -0.001 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 0.001657, max 0.033500 (between atoms 1362 and 1363)
bonds that rotated more than 30 degrees:
 atom 1 atom 2  angle  previous, current, constraint length



When I reduced this to only be the hydrogen bonds (h-bonds) I get:
---
Step 20, time 0.02 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 0.00, max 0.00 (between atoms 1322 and 1324)
bonds that rotated more than 30 degrees:
 atom 1 atom 2  angle  previous, current, constraint length
   1360   1361   37.50.1010   0.1010  0.1010
[continues with more errors until Segmentation fault 11]


I have tried both constraints with a reduced and increased deltaStep
(0.2,0.02,0.002,0.0002) no change. I have also taken constraints off, this
does work but then fails again when constraints are returned. I have tried
both Steep and cg, both fail with the same result.

The .mdp file I used (taken from the latest release of Justin¹s tutorials)
is:
; energy.mdp - used as input into grompp to generate em.tpr
; Parameters describing what to do, when to stop and what to save
integrator  = steep ; Algorithm (steep = steepest descent
minimization)
emtol   = 700.0; Stop minimization when the maximum force
< 1000.0 kJ/mol/nm
emstep  = 0.002  ; Energy step size
nsteps  = 20; Maximum number of (minimization)
steps to perform

; Parameters describing how to find the neighbors of each atom and how to
calculate the interactions
cutoff_scheme = verlet
nstlist = 20; Frequency to update the neighbor list
and long range forces
ns_type = grid  ; Method to determine neighbor list
(simple, grid)
rlist   = 1.4   ; Cut-off for making neighbor list (short
range forces)
coulombtype = PME   ; Treatment of long range electrostatic
interactions
rcoulomb= 1.4   ; Short-range electrostatic cut-off
rvdw= 1.4   ; Short-range Van der Waals cut-off
pbc = xyz   ; Periodic Boundary Conditions (yes/no)

;ADDED THIS BIT MYSELF
continuation= no   ; first dynamics run
constraint_algorithm = lincs; holonomic constraints
constraints = h-bonds; all-bonds (even heavy atom-H bonds) constrained
lincs_iter  = 1 ; accuracy of LINCS
lincs_order = 8 ; also related to accuracy


The .pdb was generated using I-TASSER. I have tried refining the hydrogen
placement via -ignh during pdb2gmx, and I also tried ³Protein Preparation²
(bond ordering - again, hydrogen bond placement) in Schrodinger. All the
same. 
 
I was hoping that I could actually watch a rough-and-ready minimisation
using Avogadro, however, the release I have keeps crashing whilst trying
to open large structures.

Any suggestions on what I could try next?


Thanks
Anthony

Dr Anthony Nash
Department of Chemistry
University College London

Skeletal Tissue Dynamics Group
Committee member of London Matrix Group @LondonMatrixGrp





On 02/10/2016 19:06, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Karel de Vries"
 wrote:

>Mark,
>
>I understand your concern.
>I need to have the hydrogens for my final answer, so an all-atom force
>field is definitely inevitable. Justification for choosing OPLS over the
>other all-atom force fields is a little more tricky, but I have found
>literature where they used OPLS-aa for similar purposes so I figured it
>should work. As for the presence of a quaternary carbon, I expect that to
>be the "alkane C", opls_139. Given that opls_135 through 138 are "alkane
>CHn", this sounds like the most reasonable interpretation of opls_139.
>
>I will try to find better documentation for the exact meaning of these
>atom
>types though, thanks.
>
>On Fri, Sep 30, 2016 at 1:49 PM, Justin Lemkul  wrote:
>
>>
>>
>> On 9/30/16 7:14 AM, Karel de Vries wrote:
>>
>>> Hi Justin,
>>>
>>> Thanks for your answer; I assumed that TopolGen is not intentionally
>>> incorrect. What I meant to ask was whether it was intentionally
>>>different
>>> from what the comments in the atomtypes.atp file indicate. I thought
>>>that
>>> perhaps, based on past experience, you know that some choices work
>>>better
>>> than others and included this into your script. I would not want to
>>>undo
>>> this by cooking up my own assignments based 

Re: [gmx-users] A charge group moved too... during backward transition of dual topologies

[Nash, Anthony]

I have tried equilibrating the non-posres simulation for an additional 2.5
ns, and taken just the final frame and velocity (so the timestep 5 5ns).
I’m also recorded a lot more to the output and I extended the fast
transition run from 100ps to 200ps so the change in softcore potentials
are not so extreme. Please note, this only happens in the reverse
transition of the dual topology, and only when I activate delta_lambda. I
end up seeing this:


Not all bonded interactions have been properly assigned to the domain
decomposition cells


A list of missing interactions:
   Angle of  14160 missing  3
 Proper Dih. of   9958 missing  3
   LJ-14 of   2214 missing  4
  exclusions of  91790 missing  8


Molecule type 'Protein2'
the first 10 missing interactions, except for exclusions:
 Proper Dih. atoms   40   43   49   52 global 63484 63487 63493
63496
   LJ-14 atoms   40   52   global 63484 63496
   Angle atoms   43   49   52  global 63487 63493 63496
 Proper Dih. atoms   44   43   49   52 global 63488 63487 63493
63496
   LJ-14 atoms   44   52   global 63488 63496
 Proper Dih. atoms   45   43   49   52 global 63489 63487 63493
63496
   LJ-14 atoms   45   52   global 63489 63496
   Angle atoms   50   49   52  global 63494 63493 63496
   Angle atoms   51   49   52  global 63495 63493 63496
   LJ-14 atoms   52   55   global 63496 63499


---
Program mdrun_mpi_d, VERSION 5.0
Source code file: 
/home/sacat2/gromacs_source/gromacs-5.0/src/gromacs/mdlib/domdec_top.c,
line: 394


Fatal error:
18 of the 130046 bonded interactions could not be calculated because some
atoms involved moved further apart than the multi-body cut-off distance
(0.934022 nm) or the two-body cut-off distance (1.42 nm), see option -rdd,
for pairs and tabulated bonds also see option -ddcheck
For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors
---



Dr Anthony Nash
Department of Chemistry
University College London





On 03/08/2016 08:11, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Nash, Anthony"
<gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
a.n...@ucl.ac.uk> wrote:

>Hi all,
>
>I¹m performing free energy calculations based on Crooks Fluctuation
>Theorem. To do this I¹ve used PMX to implement a dual topology. To keep
>things simple, it is a transmembrane polyleucine helical protein where one
>leucine is transforming into a serine - then there is the backward
>transition of serine into Leucine, required for the overlap of two
>gaussian distributions (forward and backward).
>
>I have a very complicated arrangement of steps to equilibrate the system,
>in very brief I have only included those involved once the protein is
>embedded in the bilayer (bilayer eq was around 40 ns over various stages).
>For a backward eq (serine->leucine) I used the parameters:
>
>free_energy = yes
>init_lambda = 1 ;0=forward 1=backward
>delta_lambda = 0;equilibrium
>nstdhdl = 1
>sc-coul = yes
>sc-alpha - 0.3
>sc-sigma = 0.25
>sc-power = 1
>
>
>And the eq steps:
>
>-Steepest descent (<100 kjmol threshold)
>-NVT 100ps (carbon-alpha posres 1000 kjmol)
>-NPT 100ps (carbon-alpha posres 1000 kjmol - Berendsen pressure)
>-NPT 10ns (carbon-alpha posres 1000 kjmol - PR pressure)
>-NPT 2.5 ns (turn OFF posres)
>
>The above steps produce a system with serine in the polyleu. I then
>discard the first 1.5 ns and take 10 sets of coordinates and velocities
>from the last 1ns. This is for my Œfast¹ backward (serine->leucine) 100 ps
>transformation using the MDP parameters:
>
>free_energy = yes
>init_lambda = 1 ;0=forward 1=backward
>delta_lambda = -0.2 ;transitions
>nstdhdl = 1
>sc-coul = yes
>sc-alpha - 0.3
>sc-sigma = 0.25
>sc-power = 1
>
>
>When I run the 10 Œfast¹ transition production runs approximately half of
>them immediately result in the error:
>
>Step 5:
>The charge group starting at atom 63445 moved more than the distance
>allowed by the domain decomposition (1.42) in direction Z
>distance out of cell 70.681579
>Old coordinates:  -48.465  256.312   81.655
>New coordinates:  -48.465  256.312   81.655
>Old cell boundaries in direction Z:5.728   10.974
>New cell boundaries in direction Z:5.728   10.974
>
>
>---
>Program mdrun_mpi_d, VERSION 5.0
>Source code file: 
>/home/sacat2/gromacs_source/gromacs-5.0/sr

[gmx-users] A charge group moved too... during backward transition of dual topologies

[Nash, Anthony]

Hi all,

I¹m performing free energy calculations based on Crooks Fluctuation
Theorem. To do this I¹ve used PMX to implement a dual topology. To keep
things simple, it is a transmembrane polyleucine helical protein where one
leucine is transforming into a serine - then there is the backward
transition of serine into Leucine, required for the overlap of two
gaussian distributions (forward and backward).

I have a very complicated arrangement of steps to equilibrate the system,
in very brief I have only included those involved once the protein is
embedded in the bilayer (bilayer eq was around 40 ns over various stages).
For a backward eq (serine->leucine) I used the parameters:

free_energy = yes
init_lambda = 1 ;0=forward 1=backward
delta_lambda = 0;equilibrium
nstdhdl = 1
sc-coul = yes
sc-alpha - 0.3
sc-sigma = 0.25
sc-power = 1


And the eq steps:

-Steepest descent (<100 kjmol threshold)
-NVT 100ps (carbon-alpha posres 1000 kjmol)
-NPT 100ps (carbon-alpha posres 1000 kjmol - Berendsen pressure)
-NPT 10ns (carbon-alpha posres 1000 kjmol - PR pressure)
-NPT 2.5 ns (turn OFF posres)

The above steps produce a system with serine in the polyleu. I then
discard the first 1.5 ns and take 10 sets of coordinates and velocities
from the last 1ns. This is for my Œfast¹ backward (serine->leucine) 100 ps
transformation using the MDP parameters:

free_energy = yes
init_lambda = 1 ;0=forward 1=backward
delta_lambda = -0.2 ;transitions
nstdhdl = 1
sc-coul = yes
sc-alpha - 0.3
sc-sigma = 0.25
sc-power = 1


When I run the 10 Œfast¹ transition production runs approximately half of
them immediately result in the error:

Step 5:
The charge group starting at atom 63445 moved more than the distance
allowed by the domain decomposition (1.42) in direction Z
distance out of cell 70.681579
Old coordinates:  -48.465  256.312   81.655
New coordinates:  -48.465  256.312   81.655
Old cell boundaries in direction Z:5.728   10.974
New cell boundaries in direction Z:5.728   10.974


---
Program mdrun_mpi_d, VERSION 5.0
Source code file: 
/home/sacat2/gromacs_source/gromacs-5.0/src/gromacs/mdlib/domdec.c, line:
4380


Fatal error:
A charge group moved too far between two domain decomposition steps
This usually means that your system is not well equilibrated
For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors
---


I understand what this error means, but I was wondering whether I was
doing something far more fundamentally wrong, given that every forward
transition (same eq set up, just using init_lambda = 0 and delta_lambda =
0.2) never produce this error! They always work, and monitoring my
work gives me an amazing gaussian distribution.

I have put the 2.5ns production run with no position restraints up for a
further 2.5 ns, and I¹m testing whether frames from the last ns produce
similar errors i.e., does the inclusion of serine require further
equilibration. 

I would appreciate any thoughts.

Many thanks
Anthony



-- 
Gromacs Users mailing list

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

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

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

[gmx-users] Masses and atomic radii for SASA

[Nash, Anthony]

Hi all,

I am a little concerned by the warning given (by default) when I use gmx_d
sasa Š


WARNING: Masses and atomic (Van der Waals) radii will be guessed
 based on residue and atom names, since they could not be
 definitively assigned from the information in your input
 files. These guessed numbers might deviate from the mass
 and radius of the atom type. Please check the output
 files if necessary.


NOTE: From version 5.0 gmx_d uses the Van der Waals radii
from the source below. This means the results may be different
compared to previous GROMACS versions.


 PLEASE READ AND CITE THE FOLLOWING REFERENCE 
A. Bondi
van der Waals Volumes and Radii
J. Phys. Chem. 68 (1964) pp. 441-451



How does the algorithm factor for complete bespoke atom names and atom
types in the topology? I have a post-translationally modified protein, the
particular unique area of interest does not come as standard in Gromacs so
I had to derive it from QM calculations. I now have it in my Amberffsb99
topology, but the definition is unique (atom names, atom types, residue
name etc.) Also, the warning says ³Please check the output files if
necessary.² Which ones beside the standard -o -oa and -or?

I am concerned, given the warning, that gmx_d sasa could be ignoring that
region of my protein.

Thanks
Anthony 


-- 
Gromacs Users mailing list

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

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

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

[gmx-users] Surface Distribution Function in gromacs

[Nash, Anthony]

Hi all,

A very quick sanity check question regarding one of the gromacs analysis
tools. 

Does the -surf option in gmx_d rdf (or rdf_d in 5.0.#) yield a Surface
Distribution Function plot I.e., related to the average density of water
molecules around the protein surface? I¹ve tried, and I get a plot not too
dissimilar to an RDF but with values around 0 to 300. I¹m not exactly sure
what the value on the y-xis is. There is no description in the .xvg file.

I would just like clarity that an SDF is in fact the same as this
description in the help, ³To compute the RDF with respect to the closest
position in a set in -ref instead, use -surf: if set, then -ref is
partitioned into sets based on the value of -surf, and the closest
position in each set is used."

Thanks
Anthony



-- 
Gromacs Users mailing list

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

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

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

Re: [gmx-users] minor edits to a .gro file

[Nash, Anthony]

Agreed! I have had some success using umbrella sampling, but the reaction
coordinate is quite tricky and prone to periodic boundary artefacts. I’m
putting together a script to strip out the force entries where the
corresponding frame demonstrates the long fibril protein seeing itself in
addition to its neighbour. Increasing the box, although possible, size
would make the system extremely expensive to run.

Thanks again
Anthony  



On 10/06/2016 11:25, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Mark Abraham" <gromacs.org_gmx-users-boun...@maillist.sys.kth.se
on behalf of mark.j.abra...@gmail.com> wrote:

>Hi,
>
>I certainly haven't tried something as complicated as taking three
>fragments, some charged, and morphing them to a neutral structure, but I
>would consider very seriously doing it in several stages, e.g. breaking
>into neutral fragments, dissociation of fragments, charging of fragments.
>
>Mark
>
>On Fri, Jun 10, 2016 at 12:12 PM Nash, Anthony <a.n...@ucl.ac.uk> wrote:
>
>> Hi Mark,
>>
>> I had considered (and I still am) defining a distance constraint between
>> the lysine-N and arginine-N-N on their sidechains, which defines the
>> points where the morphing into a bond formation would be.
>>
>> I am essentially trying to morph back and forth into this structure:
>> 
>>http://www.chemspider.com/Chemical-Structure.26333276.html?rid=8d8a4955-7
>>ce
>> 4-43db-b30f-12b8b6cc2f64
>> 
>><http://www.chemspider.com/Chemical-Structure.26333276.html?rid=8d8a4955-
>>7ce4-43db-b30f-12b8b6cc2f64>
>>
>>
>> Yet, I wondered whether the charges on the lysine and arginine would
>>cause
>> a bias to the result as I would be holding them very close. Still
>>thinking
>> about this one, I am going back to umbrella sampling for a while.
>>
>> Cheers
>> Anthony
>>
>>
>> Dr Anthony Nash
>> Department of Chemistry
>> University College London
>>
>>
>>
>>
>>
>> On 10/06/2016 10:37, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se
>>on
>> behalf of Mark Abraham"
>><gromacs.org_gmx-users-boun...@maillist.sys.kth.se
>> on behalf of mark.j.abra...@gmail.com> wrote:
>>
>> >Hi,
>> >
>> >One generally uses distance restraints to limit the sampling space of
>>the
>> >separated amino acids, but exactly how to implement that might take a
>>bit
>> >of iterating. GROMACS 5.1 can do inter-moleculetype restraints, but
>>that
>> >probably isn't necessary/useful in this context.
>> >
>> >Mark
>> >
>> >On Fri, Jun 10, 2016 at 11:32 AM Nash, Anthony <a.n...@ucl.ac.uk>
>>wrote:
>> >
>> >> James and Justin,
>> >>
>> >> Thanks for your suggestions. This was ultimately to hack around with
>>the
>> >> PMX toolkit for making dummy atoms in preparation for free energy
>> >> calculations. Unfortunately, I’m back to the drawing board (likely
>>to be
>> >> umbrella sampling) as an alchemical morphing of a glycated cross
>>linked
>> >> amino acid pair into two separate amino acids might work, but the
>>verse
>> >> reaction will cause all manner of constraint issues with gromacs
>>(bond
>> >> formation over a potential out of range distance).
>> >>
>> >> Thanks again
>> >> Anthony
>> >>
>> >>
>> >> Dr Anthony Nash
>> >> Department of Chemistry
>> >> University College London
>> >>
>> >>
>> >>
>> >>
>> >>
>> >> On 09/06/2016 23:52,
>>"gromacs.org_gmx-users-boun...@maillist.sys.kth.se
>> >>on
>> >> behalf of jkrie...@mrc-lmb.cam.ac.uk"
>> >> <gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
>> >> jkrie...@mrc-lmb.cam.ac.uk> wrote:
>> >>
>> >> >Plumed has a dumpatoms command (see
>> >> 
>>>http://plumed.github.io/doc-v2.1/user-doc/html/_d_u_m_p_a_t_o_m_s.html
>> >>).
>> >> >You can create virtual atoms whose position is defined based on
>> >>existing
>> >> >atoms or groups thereof. Plumed can be used as a stand-alone driver
>>or
>> >> >patched onto gromacs for on-the-fly analysis and biasing.
>> >> >
>> >> >Best wishes
>> >> >James
>> >> >
>> >> >> Hi all,
>> >> >>
>> >> >> I¹m 

[gmx-users] Missing oxyz in vmx distance output

[Nash, Anthony]

Hi all,

I¹m really hoping this is my own oversight. Using Gromacs 5.0.4 (can¹t
upgrade, this is a distribution on a cluster) I can generate each of the
output options (-oav -oall -oh -oallstat) apart from -oxyz when I run the
command:

gmx_d distance -f NPT_0_40ns_compress.gro -s NPT_0_40ns_compress.tpr -n
distance.ndx -select 'com of group "A_N" plus com of group "A_C"' -oxyz


It reports the average distance (absolute distance as there is only one
frame but I have also tried this with a .trr with the same result),
however, the distxyz.xvg is missing, even if I explicitly define ³-oxyz
distxyz.xvg².

I have also switched around the values of -pbc and -rmpbc incase there was
an undefined dependency somewhere.
 

Some of the output:
‹‹‹

GROMACS:  gmx distance, VERSION 5.0.4 (double precision)
Executable:   /usr/local/gromacs/bin/gmx_d
Library dir:  /usr/local/gromacs/share/gromacs/top
Command line:
  gmx_d distance -f NPT_0_40ns_compress.trr -s NPT_0_40ns_compress.tpr -n
distance.ndx -select 'com of group "A_N" plus com of group "A_C"' -oxyz
test.xvg


Reading file NPT_0_40ns_compress.tpr, VERSION 5.0.4 (double precision)
Reading file NPT_0_40ns_compress.tpr, VERSION 5.0.4 (double precision)
trn version: GMX_trn_file (double precision)
Last frame   2000 time 4.000
Analyzed 2001 frames, last time 4.000
A_N:
  Number of samples:  2001
  Average distance:   7.373  nm
  Standard deviation: 0.625  nm


Thanks
Anthony


Dr Anthony Nash
Department of Chemistry
University College London



-- 
Gromacs Users mailing list

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

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

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

Re: [gmx-users] minor edits to a .gro file

[Nash, Anthony]

Hi Mark,

I had considered (and I still am) defining a distance constraint between
the lysine-N and arginine-N-N on their sidechains, which defines the
points where the morphing into a bond formation would be.

I am essentially trying to morph back and forth into this structure:
http://www.chemspider.com/Chemical-Structure.26333276.html?rid=8d8a4955-7ce
4-43db-b30f-12b8b6cc2f64


Yet, I wondered whether the charges on the lysine and arginine would cause
a bias to the result as I would be holding them very close. Still thinking
about this one, I am going back to umbrella sampling for a while.

Cheers
Anthony 


Dr Anthony Nash
Department of Chemistry
University College London





On 10/06/2016 10:37, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Mark Abraham" <gromacs.org_gmx-users-boun...@maillist.sys.kth.se
on behalf of mark.j.abra...@gmail.com> wrote:

>Hi,
>
>One generally uses distance restraints to limit the sampling space of the
>separated amino acids, but exactly how to implement that might take a bit
>of iterating. GROMACS 5.1 can do inter-moleculetype restraints, but that
>probably isn't necessary/useful in this context.
>
>Mark
>
>On Fri, Jun 10, 2016 at 11:32 AM Nash, Anthony <a.n...@ucl.ac.uk> wrote:
>
>> James and Justin,
>>
>> Thanks for your suggestions. This was ultimately to hack around with the
>> PMX toolkit for making dummy atoms in preparation for free energy
>> calculations. Unfortunately, I’m back to the drawing board (likely to be
>> umbrella sampling) as an alchemical morphing of a glycated cross linked
>> amino acid pair into two separate amino acids might work, but the verse
>> reaction will cause all manner of constraint issues with gromacs (bond
>> formation over a potential out of range distance).
>>
>> Thanks again
>> Anthony
>>
>>
>> Dr Anthony Nash
>> Department of Chemistry
>> University College London
>>
>>
>>
>>
>>
>> On 09/06/2016 23:52, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se
>>on
>> behalf of jkrie...@mrc-lmb.cam.ac.uk"
>> <gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
>> jkrie...@mrc-lmb.cam.ac.uk> wrote:
>>
>> >Plumed has a dumpatoms command (see
>> >http://plumed.github.io/doc-v2.1/user-doc/html/_d_u_m_p_a_t_o_m_s.html
>>).
>> >You can create virtual atoms whose position is defined based on
>>existing
>> >atoms or groups thereof. Plumed can be used as a stand-alone driver or
>> >patched onto gromacs for on-the-fly analysis and biasing.
>> >
>> >Best wishes
>> >James
>> >
>> >> Hi all,
>> >>
>> >> I¹m looking to open a .gro file, add in six hydrogen dummy atoms (I
>>can
>> >> rename the atom name/type, I just need the correct x, y and z
>>coords) to
>> >> the end of an amino acid sidechain and save whilst preserving as
>>much of
>> >> the .gro format as it can.
>> >>
>> >>
>> >> I would normally load the crystal/derived structure as a .pdb into
>> >> Avogadro or smaller fragments from Gaussian output. Unfortunately,
>>as I
>> >> have defined a completely bespoke post-translation amino acid I can¹t
>> >> restore to .pdb with the aim of using Avogadro, too much can go
>>wrong.
>> >>
>> >>
>> >> Recommendations for Gromacs friendly editing tools would be
>>appreciated.
>> >>
>> >> Thanks
>> >> Anthony
>> >>
>> >> --
>> >> Gromacs Users mailing list
>> >>
>> >> * Please search the archive at
>> >> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>> >> posting!
>> >>
>> >> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>> >>
>> >> * For (un)subscribe requests visit
>> >> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>> >>send
>> >> a mail to gmx-users-requ...@gromacs.org.
>> >>
>> >
>>
>> --
>> Gromacs Users mailing list
>>
>> * Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>> posting!
>>
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>> * For (un)subscribe requests visit
>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>> send a mail to gmx-users-requ...@gromacs.org.
>-- 
>Gromacs Users mailing list
>
>* Please search the archive at
>http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>posting!
>
>* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
>* For (un)subscribe requests visit
>https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>send a mail to gmx-users-requ...@gromacs.org.

-- 
Gromacs Users mailing list

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

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

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

Re: [gmx-users] minor edits to a .gro file

[Nash, Anthony]

James and Justin,

Thanks for your suggestions. This was ultimately to hack around with the
PMX toolkit for making dummy atoms in preparation for free energy
calculations. Unfortunately, I’m back to the drawing board (likely to be
umbrella sampling) as an alchemical morphing of a glycated cross linked
amino acid pair into two separate amino acids might work, but the verse
reaction will cause all manner of constraint issues with gromacs (bond
formation over a potential out of range distance).

Thanks again
Anthony


Dr Anthony Nash
Department of Chemistry
University College London





On 09/06/2016 23:52, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of jkrie...@mrc-lmb.cam.ac.uk"
 wrote:

>Plumed has a dumpatoms command (see
>http://plumed.github.io/doc-v2.1/user-doc/html/_d_u_m_p_a_t_o_m_s.html ).
>You can create virtual atoms whose position is defined based on existing
>atoms or groups thereof. Plumed can be used as a stand-alone driver or
>patched onto gromacs for on-the-fly analysis and biasing.
>
>Best wishes
>James
>
>> Hi all,
>>
>> I¹m looking to open a .gro file, add in six hydrogen dummy atoms (I can
>> rename the atom name/type, I just need the correct x, y and z coords) to
>> the end of an amino acid sidechain and save whilst preserving as much of
>> the .gro format as it can.
>>
>>
>> I would normally load the crystal/derived structure as a .pdb into
>> Avogadro or smaller fragments from Gaussian output. Unfortunately, as I
>> have defined a completely bespoke post-translation amino acid I can¹t
>> restore to .pdb with the aim of using Avogadro, too much can go wrong.
>>
>>
>> Recommendations for Gromacs friendly editing tools would be appreciated.
>>
>> Thanks
>> Anthony
>>
>> --
>> Gromacs Users mailing list
>>
>> * Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>> posting!
>>
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>> * For (un)subscribe requests visit
>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>>send
>> a mail to gmx-users-requ...@gromacs.org.
>>
>

-- 
Gromacs Users mailing list

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

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

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

[gmx-users] minor edits to a .gro file

[Nash, Anthony]

Hi all,

I¹m looking to open a .gro file, add in six hydrogen dummy atoms (I can
rename the atom name/type, I just need the correct x, y and z coords) to
the end of an amino acid sidechain and save whilst preserving as much of
the .gro format as it can.


I would normally load the crystal/derived structure as a .pdb into
Avogadro or smaller fragments from Gaussian output. Unfortunately, as I
have defined a completely bespoke post-translation amino acid I can¹t
restore to .pdb with the aim of using Avogadro, too much can go wrong.
 

Recommendations for Gromacs friendly editing tools would be appreciated.

Thanks
Anthony 

-- 
Gromacs Users mailing list

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

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

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

Re: [gmx-users] Clarity on TI free energy terms

[Nash, Anthony]

Dear Hannes,

Thanks for all the help yesterday, it helped. I, hopefully, have just the
one final question.

I am still a little confused how Gromacs deals with the interactions (vdW
& Coul) with the environment when a soft-core potential has been used to
switch between molecules (typeA and typeB in the topology file). I
understand how lambda can be used to 1) phase out the charge then 2) phase
out the vdW interactions for a molecule that is disappearing (or reverse
for appearing) to capture hydration energetics e.g., Justin¹s methane in
water FE tutorial. 

However, in the case of a D2K (aspartic-acid to lysine) dual-topology, how
are the vdW and charges brought back into play for typeB, when the
interactions of typeA have been coupled to lambda as per your example
below? 

fep-lambdas = 0.0 0.2 0.4 0.6 0.8 1.0  1.0  1.0  etc.
vdw-lambdas = 0.0 0.0 0.0 0.0 0.0 0.0  0.05 0.10 etc.
mass-lambdas as above

I apologise if I may have misunderstood one of your earlier explanations,
I think this is the only piece of the puzzle left for me to understand.

Many thanks
Anthony



On 01/06/2016 19:04, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Hannes Loeffler"
<gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
hannes.loeff...@stfc.ac.uk> wrote:

>On Wed, 1 Jun 2016 16:15:31 +
>"Nash, Anthony" <a.n...@ucl.ac.uk> wrote:
>
>> In the mean while, do you know of any tutorials (besides the methane
>> in water FE tutorial) regarding TI for amino acid substitution?
>
>I am not aware of one.  You could try
>http://www.alchemistry.org/
>
>
>> And by ³q_off² and ³vdw_on/off², I assume you are referring to the
>> Œvalue¹ of the couple-lambda0: parameter?
>
>No.  I'm referring to the respective lambda paths.
>
>
>Cheers,
>Hannes.
>-- 
>Gromacs Users mailing list
>
>* Please search the archive at
>http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>posting!
>
>* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
>* For (un)subscribe requests visit
>https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>send a mail to gmx-users-requ...@gromacs.org.

-- 
Gromacs Users mailing list

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

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

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

Re: [gmx-users] Clarity on TI free energy terms

[Nash, Anthony]

In the mean while, do you know of any tutorials (besides the methane in
water FE tutorial) regarding TI for amino acid substitution? And by
“q_off” and “vdw_on/off”, I assume you are referring to the ‘value’ of the
couple-lambda0: parameter?

Thanks
Anthony

Dr Anthony Nash
Department of Chemistry
University College London





On 01/06/2016 16:55, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Nash, Anthony"
<gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
a.n...@ucl.ac.uk> wrote:

>Thanks for all of this material. I’ll take some time and digest what
>you’ve said. 
>
>No doubt I’ll have a few more questions tomorrow ;-)
>
>Thanks
>Anthony 
>
>
>On 01/06/2016 16:38, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
>behalf of Hannes Loeffler"
><gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
>hannes.loeff...@stfc.ac.uk> wrote:
>
>>On Wed, 1 Jun 2016 15:00:51 +
>>"Nash, Anthony" <a.n...@ucl.ac.uk> wrote:
>>
>>> >  This also assumes that
>>> >you have vanishing atoms only.  If you have appearing atoms only you
>>> >would obviously have to revers the order, and when you have both you
>>> >will have to run with two mdp/tpr setups.
>>> 
>>> 
>>> With aspartic acid transforming into lysine on one polypeptide chain
>>> and then lysine transforming into aspartic acid polypeptide chain,
>>> all in the same system (keeps the charges the same and is a complete
>>> thermodynamic cycle), I will have both appearing and disappearing
>>> atoms. How do you mean ³to run with two mdp/tpr setups²? Is there an
>>> example you could give (which I am grateful for) or one in the manual?
>>
>>Lambda paths are not selective in the sense that you could apply them
>>to only a subset of the system.  So if you have both disappearing and
>>appearing atoms you have to:
>>1) turn off the charges on the disappearing group (or alternatively all
>>charges to avoid a charged total system), q_off
>>2) turn off the vdW parameters for the disappearing group, vdw_off
>>3) turn on the vdW parameters for the appearing group, vdw_on
>>4) turn on the charges of the appearing group (or all charges), q_on
>>
>>If you try to do this with Gromacs you will realise that the best you
>>can do is combine this into 2 mdp steps: 1) q_off and vdw_on/off and 2)
>>q_on.  Alternatively you can combine the vdw bit with 2) but that
>>doesn't make any difference.  Of course, if you could assume that a
>>combined Coulomb/vdW transformation will work, this would be a
>>non-issue...
>>
>>Also, I wonder how pmx maps aspartate onto lysine.  I would think that
>>the gamma-carboxylate should not map onto the gamma-methylene but
>>rather the residue should be branched off at the C-beta and so
>>duplicate the rest of the residue.
>>-- 
>>Gromacs Users mailing list
>>
>>* Please search the archive at
>>http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>>posting!
>>
>>* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>>* For (un)subscribe requests visit
>>https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>>send a mail to gmx-users-requ...@gromacs.org.
>
>-- 
>Gromacs Users mailing list
>
>* Please search the archive at
>http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>posting!
>
>* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
>* For (un)subscribe requests visit
>https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>send a mail to gmx-users-requ...@gromacs.org.

-- 
Gromacs Users mailing list

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

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

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

Re: [gmx-users] Clarity on TI free energy terms

[Nash, Anthony]

Thanks for all of this material. I’ll take some time and digest what
you’ve said. 

No doubt I’ll have a few more questions tomorrow ;-)

Thanks
Anthony 


On 01/06/2016 16:38, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Hannes Loeffler"
<gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
hannes.loeff...@stfc.ac.uk> wrote:

>On Wed, 1 Jun 2016 15:00:51 +0000
>"Nash, Anthony" <a.n...@ucl.ac.uk> wrote:
>
>> >  This also assumes that
>> >you have vanishing atoms only.  If you have appearing atoms only you
>> >would obviously have to revers the order, and when you have both you
>> >will have to run with two mdp/tpr setups.
>> 
>> 
>> With aspartic acid transforming into lysine on one polypeptide chain
>> and then lysine transforming into aspartic acid polypeptide chain,
>> all in the same system (keeps the charges the same and is a complete
>> thermodynamic cycle), I will have both appearing and disappearing
>> atoms. How do you mean ³to run with two mdp/tpr setups²? Is there an
>> example you could give (which I am grateful for) or one in the manual?
>
>Lambda paths are not selective in the sense that you could apply them
>to only a subset of the system.  So if you have both disappearing and
>appearing atoms you have to:
>1) turn off the charges on the disappearing group (or alternatively all
>charges to avoid a charged total system), q_off
>2) turn off the vdW parameters for the disappearing group, vdw_off
>3) turn on the vdW parameters for the appearing group, vdw_on
>4) turn on the charges of the appearing group (or all charges), q_on
>
>If you try to do this with Gromacs you will realise that the best you
>can do is combine this into 2 mdp steps: 1) q_off and vdw_on/off and 2)
>q_on.  Alternatively you can combine the vdw bit with 2) but that
>doesn't make any difference.  Of course, if you could assume that a
>combined Coulomb/vdW transformation will work, this would be a
>non-issue...
>
>Also, I wonder how pmx maps aspartate onto lysine.  I would think that
>the gamma-carboxylate should not map onto the gamma-methylene but
>rather the residue should be branched off at the C-beta and so
>duplicate the rest of the residue.
>-- 
>Gromacs Users mailing list
>
>* Please search the archive at
>http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>posting!
>
>* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
>* For (un)subscribe requests visit
>https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>send a mail to gmx-users-requ...@gromacs.org.

-- 
Gromacs Users mailing list

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

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

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

Re: [gmx-users] Clarity on TI free energy terms

[Nash, Anthony]

Dear Hannes, please see my comment below..


On 01/06/2016 14:45, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Hannes Loeffler"
<gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
hannes.loeff...@stfc.ac.uk> wrote:

>On Wed, 1 Jun 2016 12:06:20 +0000
>"Nash, Anthony" <a.n...@ucl.ac.uk> wrote:
>
>> vdw_lambdas  = 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35
>> 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00
>> mass_lambdas = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
>> 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
>> fep_lambdas =  0.0 0.1 0.1 0.2 0.2 0.3 0.3 0.4 0.4 0.5 0.5 0.6 0.6
>> 0.7 0.7 0.8 0.8 0.9 0.9 1.0 1.0
>
>This will transform both vdW _and_ Coulomb terms at the same time but
>at a different pace.  Maybe you want something like
>
>fep-lambdas = 0.0 0.2 0.4 0.6 0.8 1.0  1.0  1.0  etc.
>vdw-lambdas = 0.0 0.0 0.0 0.0 0.0 0.0  0.05 0.10 etc.
>mass-lambdas as above
>
>This will first transform Coulomb and bonded terms in the first six
>lambdas and vdW from the 7th lambda onwards.


That is beginning to make perfect sense now. Many thanks for that help.



>  This also assumes that
>you have vanishing atoms only.  If you have appearing atoms only you
>would obviously have to revers the order, and when you have both you
>will have to run with two mdp/tpr setups.


With aspartic acid transforming into lysine on one polypeptide chain and
then lysine transforming into aspartic acid polypeptide chain, all in the
same system (keeps the charges the same and is a complete thermodynamic
cycle), I will have both appearing and disappearing atoms. How do you mean
³to run with two mdp/tpr setups²? Is there an example you could give
(which I am grateful for) or one in the manual?


>> couple-moltype   = protein
>> couple-lambda0   = vdw-q
>> couple-lambda1   = vdw-q
>
>This will decouple all atoms in the 'protein' selection from the
>environment.  This is for absolute transformation and not what you seem
>to want to do i.e. a relative transformation of one residue into
>another.  So avoid those parameters.

Thanks for that information.


Many thanks
Anthony

-- 
Gromacs Users mailing list

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

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

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

Re: [gmx-users] Clarity on TI free energy terms

[Nash, Anthony]

Hi Hannes,

Thanks for the reply. At the moment for a change in single amino acid in a
triplet (a pair of triplets, showing forward and reverse change) I am
settling with:

vdw_lambdas  = 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40
0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00
mass_lambdas = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0
fep_lambdas =  0.0 0.1 0.1 0.2 0.2 0.3 0.3 0.4 0.4 0.5 0.5 0.6 0.6 0.7 0.7
0.8 0.8 0.9 0.9 1.0 1.0


free_energy  = yes
init_lambda_state= 0
delta_lambda = 0
calc_lambda_neighbors= 1; only immediate neighboring windows


; Options for the decoupling
sc-alpha = 0.5
sc-coul  = yes
sc-power = 1.0
sc-sigma = 0.3
couple-moltype   = protein
couple-lambda0   = vdw-q
couple-lambda1   = vdw-q
couple-intramol  = no
nstdhdl  = 10



The  mass will be conserved (it is a full cycle in one system). Everything
else scales with feb_lambdas apart from vdw which will (I am guessing)
require more sampling. This is my first attempt, I don't expect to get a
true understanding of it yet :-)

Thanks again
Anthony





On 01/06/2016 12:43, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Hannes Loeffler"
<gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
hannes.loeff...@stfc.ac.uk> wrote:

>
>Set the vector to all-zeroes (or ones).
>
>
>On Wed, 1 Jun 2016 09:47:59 +
>"Nash, Anthony" <a.n...@ucl.ac.uk> wrote:
>
>> Hi Hannes,
>> 
>> 
>> Many thanks for the reply. With regards to your final comment I
>> understand conserving mass in theory, but I am a little confused
>> regarding, "keep the mass-lambdas at one end-point as they can
>> interact badly with constraints". I am testing pmx on a two-molecule
>> one-system I.e., G-D2K-G and G-K2D-G in the same system. How ought I
>> define the mass-lambdas for this system? (nothing accurate, just an
>> example would be great)?
>> 
>> Thanks
>> Anthony
>> 
>> 
>> On 01/06/2016 09:55,
>> "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
>> Hannes Loeffler" <gromacs.org_gmx-users-boun...@maillist.sys.kth.se
>> on behalf of hannes.loeff...@stfc.ac.uk> wrote:
>> 
>> >On Wed, 1 Jun 2016 07:54:56 +
>> >"Nash, Anthony" <a.n...@ucl.ac.uk> wrote:
>> >
>> >> In the tutorial, charges are off in the topology and the
>> >> electrostatic coupling to lambda remains 0 throughout the 20
>> >> windows. I assume setting col_lambdas=0 0 0 Š was for that very
>> >> reason I.e., the charges were off? Could the charges not have been
>> >> left on and col_lambdas defined similar to vdw_lambdas?
>> >> (I understand that if charges remain constant, as vdw turns off,
>> >> the system will probably blow up as attraction brings molecules
>> >> infinitely close).
>> >
>> >Technically, Gromacs allows you to vary both vdW and Coulomb lambdas
>> >simultaneously because Gromacs can apply softcore potentials to both.
>> >In practice though it seems that many workers still prefer to
>> >separate the two terms from each other.
>> >
>> > 
>> >> If my transition is from a small molecule into a small molecule
>> >> e.g., G-D-G to G-K-D, (the PMX paper) should I define all three
>> >> lambdas: vdw_lambdas, col_lambdas and bonds_lambdas? Between
>> >> states A and B, VdW, charges and bonds are all changing.
>> >
>> >Lambda paths are only about separating the various force field terms
>> >from each other.  If you do not explicitly state any of those lambda
>> >vectors they will adopt they same lambdas as specified in
>> >fep-lambdas, see manual.  I do not see a reason why you would want
>> >to separate out the bonded terms as well.  They are subject to a
>> >linear transformation only anyway.
>> >
>> >What you may want to do is to keep the mass-lambdas at one end-point
>> >as they can interact badly with constraints.  In a proper
>> >thermodynamic cycle mass contributions must perfectly cancel.
>> >
>> >
>> >Cheers,
>> >Hannes.
>> >-- 
>> >Gromacs Users mailing list
>> >
>> >* Please search the archive at
>> >http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>> >posting!
>> >
>> >* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>> >
>> &g

Re: [gmx-users] Clarity on TI free energy terms

[Nash, Anthony]

Hi Hannes,


Many thanks for the reply. With regards to your final comment I understand
conserving mass in theory, but I am a little confused regarding, "keep the
mass-lambdas at one end-point as they can interact badly with
constraints". I am testing pmx on a two-molecule one-system I.e., G-D2K-G
and G-K2D-G in the same system. How ought I define the mass-lambdas for
this system? (nothing accurate, just an example would be great)?

Thanks
Anthony


On 01/06/2016 09:55, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Hannes Loeffler"
<gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
hannes.loeff...@stfc.ac.uk> wrote:

>On Wed, 1 Jun 2016 07:54:56 +
>"Nash, Anthony" <a.n...@ucl.ac.uk> wrote:
>
>> In the tutorial, charges are off in the topology and the electrostatic
>> coupling to lambda remains 0 throughout the 20 windows. I assume
>> setting col_lambdas=0 0 0 Š was for that very reason I.e., the
>> charges were off? Could the charges not have been left on and
>> col_lambdas defined similar to vdw_lambdas?
>> (I understand that if charges remain constant, as vdw turns off, the
>> system will probably blow up as attraction brings molecules infinitely
>> close).
>
>Technically, Gromacs allows you to vary both vdW and Coulomb lambdas
>simultaneously because Gromacs can apply softcore potentials to both.
>In practice though it seems that many workers still prefer to separate
>the two terms from each other.
>
> 
>> If my transition is from a small molecule into a small molecule e.g.,
>> G-D-G to G-K-D, (the PMX paper) should I define all three lambdas:
>> vdw_lambdas, col_lambdas and bonds_lambdas? Between states A and B,
>> VdW, charges and bonds are all changing.
>
>Lambda paths are only about separating the various force field terms
>from each other.  If you do not explicitly state any of those lambda
>vectors they will adopt they same lambdas as specified in fep-lambdas,
>see manual.  I do not see a reason why you would want to separate out
>the bonded terms as well.  They are subject to a linear transformation
>only anyway.
>
>What you may want to do is to keep the mass-lambdas at one end-point as
>they can interact badly with constraints.  In a proper thermodynamic
>cycle mass contributions must perfectly cancel.
>
>
>Cheers,
>Hannes.
>-- 
>Gromacs Users mailing list
>
>* Please search the archive at
>http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>posting!
>
>* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
>* For (un)subscribe requests visit
>https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>send a mail to gmx-users-requ...@gromacs.org.

-- 
Gromacs Users mailing list

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

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

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

[gmx-users] Clarity on TI free energy terms

[Nash, Anthony]

Dear all,

I¹m trying to understand the finesse behind the TI free energy in gromacs,
before taking it anywhere near a real production run, by running through
the FE methane in solvent tutorial and the thermodynamic cycles of small
peptides in the PMX paper. I roughly-understand a fair chunk, however, I
would appreciate some clarity on one or two points.

In the tutorial, charges are off in the topology and the electrostatic
coupling to lambda remains 0 throughout the 20 windows. I assume setting
col_lambdas=0 0 0 Š was for that very reason I.e., the charges were off?
Could the charges not have been left on and col_lambdas defined similar to
vdw_lambdas? 
(I understand that if charges remain constant, as vdw turns off, the
system will probably blow up as attraction brings molecules infinitely
close). 

If my transition is from a small molecule into a small molecule e.g.,
G-D-G to G-K-D, (the PMX paper) should I define all three lambdas:
vdw_lambdas, col_lambdas and bonds_lambdas? Between states A and B, VdW,
charges and bonds are all changing.

Many thanks
Anthony

-- 
Gromacs Users mailing list

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

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

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

Re: [gmx-users] Free Energy Topology between A and B.

[Nash, Anthony]

So there is! Many thanks for bringing this to my attention.


Thanks
Anthony 

On 29/05/2016 20:40, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Justin Lemkul"
<gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
jalem...@vt.edu> wrote:

>
>
>On 5/29/16 3:37 PM, Nash, Anthony wrote:
>> Dear all,
>>
>> I¹m a total newbie when it comes to Thermodynamic Integration, and until
>> now I¹ve been happy with umbrella sampling. However, I¹ve found myself
>>in
>> a situation where I believe TI would be the most appropriate technique.
>>
>> I would like to determine the energetic contribution that a mutant amino
>> acid has to a protein-protein interaction. The structure is very
>>complex,
>> so pulling along a reaction pathway is out of the question, thus an
>> alchemical path might work.
>>
>> I wish to integrate from state A to state B. The number of atoms in the
>> wild type amino acid of state A are fewer that those in the mutant amino
>> acid of state B. The example given in section 5.7.4 of the manual
>> demonstrate integrating from propanol to pentane, but due to the
>>GROMOS-96
>> FF each topology shares the same number of atoms. Atom PC2 and PC3 of
>> state A does not require redefining in typeB chargeB and massB, they are
>> transferred across as they are in both structures. How does one define a
>> mixed A B topology when the number of atoms differ? I.e., if those two
>> atoms in the above example were not needed in state B.
>>
>
>You need to define transformations between dummy and real atoms.  Look
>into the 
>pmx program (this was mentioned on the list just the other day for an
>identical 
>problem).
>
>-Justin
>
>-- 
>==
>
>Justin A. Lemkul, Ph.D.
>Ruth L. Kirschstein NRSA Postdoctoral Fellow
>
>Department of Pharmaceutical Sciences
>School of Pharmacy
>Health Sciences Facility II, Room 629
>University of Maryland, Baltimore
>20 Penn St.
>Baltimore, MD 21201
>
>jalem...@outerbanks.umaryland.edu | (410) 706-7441
>http://mackerell.umaryland.edu/~jalemkul
>
>==
>-- 
>Gromacs Users mailing list
>
>* Please search the archive at
>http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>posting!
>
>* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
>* For (un)subscribe requests visit
>https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>send a mail to gmx-users-requ...@gromacs.org.

-- 
Gromacs Users mailing list

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

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

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

[gmx-users] Free Energy Topology between A and B.

[Nash, Anthony]

Dear all,

I¹m a total newbie when it comes to Thermodynamic Integration, and until
now I¹ve been happy with umbrella sampling. However, I¹ve found myself in
a situation where I believe TI would be the most appropriate technique.

I would like to determine the energetic contribution that a mutant amino
acid has to a protein-protein interaction. The structure is very complex,
so pulling along a reaction pathway is out of the question, thus an
alchemical path might work.

I wish to integrate from state A to state B. The number of atoms in the
wild type amino acid of state A are fewer that those in the mutant amino
acid of state B. The example given in section 5.7.4 of the manual
demonstrate integrating from propanol to pentane, but due to the GROMOS-96
FF each topology shares the same number of atoms. Atom PC2 and PC3 of
state A does not require redefining in typeB chargeB and massB, they are
transferred across as they are in both structures. How does one define a
mixed A B topology when the number of atoms differ? I.e., if those two
atoms in the above example were not needed in state B.

Thanks
Anthony

-- 
Gromacs Users mailing list

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

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

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

Re: [gmx-users] gmx hbond - specify precise atom names involved

[Nash, Anthony]

Thanks Justin,


I’ll give that a try. I assume this approach would still require the .trr
to be converted to a .gro file, and then the customised names
‘search-replaced’ with the respective ‘fake’ names?

Thanks
Anthony

On 03/05/2016 16:22, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Justin Lemkul"
<gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
jalem...@vt.edu> wrote:

>
>
>On 5/3/16 9:16 AM, Nash, Anthony wrote:
>>
>> Hi all,
>>
>> Can gmx hbond accept user specified atoms for the donors (default OH and
>> NH) and acceptor (default O and N)? I don¹t seem to find any mention of
>> this in the -help text.
>>
>
>It's hard-coded in the source, so it's rather inflexible.
>
>> I have a post-trans modified protein from a rather bulk cross-linked
>> peptide chain. I defined unique atom times but I have used a unique set
>>of
>> atom names.
>>
>
>Create a "fake" topology that uses O,N,H as the first character in
>applicable 
>atom names.  Then run the analysis using this .tpr file.  We did this a
>while 
>back with, e.g. thiols and a few others.
>
>-Justin
>
>-- 
>==
>
>Justin A. Lemkul, Ph.D.
>Ruth L. Kirschstein NRSA Postdoctoral Fellow
>
>Department of Pharmaceutical Sciences
>School of Pharmacy
>Health Sciences Facility II, Room 629
>University of Maryland, Baltimore
>20 Penn St.
>Baltimore, MD 21201
>
>jalem...@outerbanks.umaryland.edu | (410) 706-7441
>http://mackerell.umaryland.edu/~jalemkul
>
>==
>-- 
>Gromacs Users mailing list
>
>* Please search the archive at
>http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>posting!
>
>* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
>* For (un)subscribe requests visit
>https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>send a mail to gmx-users-requ...@gromacs.org.

-- 
Gromacs Users mailing list

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

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

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

[gmx-users] gmx hbond - specify precise atom names involved

[Nash, Anthony]

Hi all,

Can gmx hbond accept user specified atoms for the donors (default OH and
NH) and acceptor (default O and N)? I don¹t seem to find any mention of
this in the -help text.

I have a post-trans modified protein from a rather bulk cross-linked
peptide chain. I defined unique atom times but I have used a unique set of
atom names. 

Thanks
Anthony 




>

-- 
Gromacs Users mailing list

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

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

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

[gmx-users] Constant Density

[Nash, Anthony]

Hi all,

At the risk of bending the rules of thermodynamics, I¹m wondering whether
Gromacs can maintain density of a water box (0.750 g/L density of water in
a collagen fibril environment) whilst applying an NPT ensemble?

gmx_d solvate, fills up to 2/3 of my truncated oct cell, with my protein
at the centre. An NVT simulation does not redistribute the water. To do
this I need to perform an NPT run, but even so this rapidly shrinks the
box. 

I have considering looking up the necessary pressure value for this
density, however I¹m a bit uncertain whether this maintains any
physiological realism for the protein itself.

I¹ve tried googling with little success.

Any thoughts are appreciated.
Thanks
Anthony 

-- 
Gromacs Users mailing list

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

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

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

Re: [gmx-users] Thermodynamic integration

[Nash, Anthony]

Hi Mark,

I will give that a thorough read. I was wondering if you could possibly
comment on whether TI is an appropriate tool for calculating the free
energy difference between two states, A―> non-glycated side chain, and b―>
glycated side chain? Most examples given focus on the inclusion/exclusion
of some small molecule in some large system.

I have tried umbrella sampling and although my results are extremely
interesting, I’ve had to manipulate the initial systems to take it from a
crystal structure with defined periodicity in x-y-z dimensions, to a slab
in the x-y plane, and a water bath in the z plane.
 

Thanks
Anthony 




On 18/04/2016 11:53, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Mark Abraham" <gromacs.org_gmx-users-boun...@maillist.sys.kth.se
on behalf of mark.j.abra...@gmail.com> wrote:

>Hi,
>
>Also you might consider pmx
>http://www.ncbi.nlm.nih.gov/pmc/articles/PMC4365728/ for such topology
>generation. There is further work in the pipeline, so do get in touch with
>Bert if there's something of interest.
>
>Mark
>
>On Mon, Apr 18, 2016 at 11:56 AM Nash, Anthony <a.n...@ucl.ac.uk> wrote:
>
>> From the site, “..or the free energy of a mutation of a side chain.”
>>
>> I think this is what I am after. Many thanks for the link.
>>
>> Anthony
>>
>>
>>
>> On 18/04/2016 10:42, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se
>>on
>> behalf of Hannes Loeffler"
>> <gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
>> hannes.loeff...@stfc.ac.uk> wrote:
>>
>> >A good starting point is http://www.alchemistry.org/ which has quite a
>> >lot of detail on relative alchemical free energy simulations (not only
>> >TI).
>> >
>> >On Mon, 18 Apr 2016 09:27:02 +
>> >"Nash, Anthony" <a.n...@ucl.ac.uk> wrote:
>> >
>> >> Hi all,
>> >>
>> >> I¹m looking for a guide on performing TI between a protein in its
>> >> crystal periodicity with a particular residue (state A), to the same
>> >> system but with a different residue (state B).
>> >>
>> >> I¹m currently using
>> >>
>> >>
>> 
>>http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/fr
>> >>ee
>> >> _energy/01_theory.html as a guide, however this is more, if I
>> >> understand having read through it, on the presence-to-absence of
>> >> methane in a water solvent rather than a replacement with something
>> >> else.
>> >>
>> >> I haven¹t had too much luck googling and I¹m looking piecemeal
>> >> through the manual with little success.
>> >>
>> >> Thanks
>> >> Anthony
>> >>
>> >>
>> >>
>> >
>> >--
>> >Gromacs Users mailing list
>> >
>> >* Please search the archive at
>> >http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>> >posting!
>> >
>> >* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>> >
>> >* For (un)subscribe requests visit
>> >https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>> >send a mail to gmx-users-requ...@gromacs.org.
>>
>> --
>> Gromacs Users mailing list
>>
>> * Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>> posting!
>>
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>> * For (un)subscribe requests visit
>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>> send a mail to gmx-users-requ...@gromacs.org.
>-- 
>Gromacs Users mailing list
>
>* Please search the archive at
>http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>posting!
>
>* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
>* For (un)subscribe requests visit
>https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>send a mail to gmx-users-requ...@gromacs.org.

-- 
Gromacs Users mailing list

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

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

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

[gmx-users] Thermodynamic integration

[Nash, Anthony]

Hi all,

I¹m looking for a guide on performing TI between a protein in its crystal
periodicity with a particular residue (state A), to the same system but
with a different residue (state B).

I¹m currently using
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/free
_energy/01_theory.html as a guide, however this is more, if I understand
having read through it, on the presence-to-absence of methane in a water
solvent rather than a replacement with something else.

I haven¹t had too much luck googling and I¹m looking piecemeal through the
manual with little success.

Thanks
Anthony  



-- 
Gromacs Users mailing list

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

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

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

[gmx-users] pull code for Gromacs 5

[Nash, Anthony]

Hi all,

I¹m very unfamiliar with the pull code as of Gromacs 5. Unfortunately my
system is not experiencing any noticeable Œpull¹. From the options below,
which is the group experiencing the pull and which is the reference group?
 Would applying a set of position restraints on the reference group
prevent (in theory) the pulling group to move? From someone only having
ran pull code on earlier versions of gromacs, am I missing anything
blindingly obvious?


pull= umbrella
pull-geometry   = direction-periodic
pull-dim= N N Y
pull-start  = yes
pull-ncoords= 1
pull-ngroups= 2
pull-group1-name= fix_collagen
pull-group1-pbcatom = 0
pull-group2-name= pull_group
pull-coord1-groups  = 1 2
pull-coord1-rate= 0.1
pull-coord1-k   = 1000



Context: I am pulling an explicit collagen molecule away from its
neighbour. They are in a very tight Œrectangle¹, thus they are
experiencing a pseudo fibril environment. The need for direction-periodic
is due to the length I am pulling the molecule (past half the box size).

Thanks
Anthony



Dr Anthony Nash
Department of Chemistry
University College London



-- 
Gromacs Users mailing list

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

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

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

Re: [gmx-users] Can I fix two of the cubic cell dimensions during an NPT simulation?

[Nash, Anthony]

Justin, that’s awesome. Thanks



On 14/03/2016 22:34, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Justin Lemkul"
<gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
jalem...@vt.edu> wrote:

>
>
>On 3/14/16 6:31 PM, Nash, Anthony wrote:
>> Hi all,
>>
>> Is there a way of keeping the x, y box dimensions fixed during an NPT
>> simulation, with changes to volume only changing in the Z dimension?
>> Semiisotropic is not quite working out, see below.
>>
>>
>>
>> Context: I want a coiled-coil dimer aligned in the Z direction. Each
>> coiled-coil will see it¹s explicit neighbour, plus all period images.
>>This
>> will give an impression of a collagen fibrillar environment. Each
>>protein
>> should not be be able to tilt off the z-axis, due to the packed nature.
>> Either end of the two coiled coils will be Œbulk¹ water which can expand
>> and contract from changes in volume. If the x or y dimension changes
>>there
>> is a risk to water forming a layer between an explicit coiled-coil and
>>one
>> of its period images, thus losing the fibrillar environment effect.
>>
>
>Set compressibility = 0 for x-y.
>
>-Justin
>
>-- 
>==
>
>Justin A. Lemkul, Ph.D.
>Ruth L. Kirschstein NRSA Postdoctoral Fellow
>
>Department of Pharmaceutical Sciences
>School of Pharmacy
>Health Sciences Facility II, Room 629
>University of Maryland, Baltimore
>20 Penn St.
>Baltimore, MD 21201
>
>jalem...@outerbanks.umaryland.edu | (410) 706-7441
>http://mackerell.umaryland.edu/~jalemkul
>
>==
>-- 
>Gromacs Users mailing list
>
>* Please search the archive at
>http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>posting!
>
>* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
>* For (un)subscribe requests visit
>https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>send a mail to gmx-users-requ...@gromacs.org.

-- 
Gromacs Users mailing list

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

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

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

[gmx-users] Can I fix two of the cubic cell dimensions during an NPT simulation?

[Nash, Anthony]

Hi all, 

Is there a way of keeping the x, y box dimensions fixed during an NPT
simulation, with changes to volume only changing in the Z dimension?
Semiisotropic is not quite working out, see below.



Context: I want a coiled-coil dimer aligned in the Z direction. Each
coiled-coil will see it¹s explicit neighbour, plus all period images. This
will give an impression of a collagen fibrillar environment. Each protein
should not be be able to tilt off the z-axis, due to the packed nature.
Either end of the two coiled coils will be Œbulk¹ water which can expand
and contract from changes in volume. If the x or y dimension changes there
is a risk to water forming a layer between an explicit coiled-coil and one
of its period images, thus losing the fibrillar environment effect.

Thanks
Anthony 

Dr Anthony Nash
Department of Chemistry
University College London

-- 
Gromacs Users mailing list

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

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

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

[gmx-users] Suggestions on running simulations of very long polypeptide chains

[Nash, Anthony]

Hi all,

I¹m looking to run MD simulations of regions of a collagen molecule. A
whole collegen molecule is made up of three polypeptide chains, each
around 1000 residues long (gross generalisation as there are around 24
different collagen protein families). I am only interested in modelling a
section of 30 residues in length, and then as a dimer (30 per chain, three
chains per molecule, 2 molecules to make the dimer, 180 residues in the
whole system)

I tried looking at the structure as a free-floating dimer in solvent. Two
things occurred, 1) they became monomeric, 2) the system box had to be
over twice the size in all three dimensions to consider all of the
possible conformations without seeing its period image. This went from
being a very small system to a very large system due to the solvent.

After some careful consideration I think my best approach would be to turn
each 30 long residue into a polypeptide chain which covalently binds back
onto itself, giving the illusion of a very long polypeptide chain.  Are
there any tutorials of best practice on this techniques, and in particular
how to manage the pressure on the adjusting volume?

Many thanks
Anthony

Dr Anthony Nash
Department of Chemistry
University College London

-- 
Gromacs Users mailing list

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

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

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

Re: [gmx-users] Minimising forces for vibrational normal mode analysis

[Nash, Anthony]

Thanks Peter, I thought as much. At the time it just seemed like an odd
bit of information to report.

Many thanks
Anthony 

Dr Anthony Nash
Department of Chemistry
University College London





On 01/03/2016 03:15, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Peter Stern" <gromacs.org_gmx-users-boun...@maillist.sys.kth.se
on behalf of peter.st...@weizmann.ac.il> wrote:

>This is a non-problem.  You will never get zero eigenvalues for such a
>large system since you are not at an absolute minimum.  The fact that you
>got no negative eigenvalues and the six lowest eigenvalues are very small
>(surely much, much smaller than the subsequent eigenvalues) indicates
>that you are at a reasonably good minimum. But don't be misled. Another
>(different) minimization could bring the system to a different minimum
>(also good) with somewhat different eigenvalues and eigenvectors.
>
>Regards,
>Peter Stern
>
>Sent from my iPhone
>
>> On 29 Feb 2016, at 3:23 PM, Nash, Anthony <a.n...@ucl.ac.uk> wrote:
>> 
>> Hi all,
>> 
>> Hoping this will be the final stumbling block. I’ve successfully
>>minimised
>> my structure to a a very small max force:
>> Potential Energy  = -1.93736854460113e+03
>> Maximum force =  7.39399083846082e-04 on atom 34
>> Norm of force =  1.64068483698544e-04
>> 
>> 
>> When I perform integrator=nm there are no warning of potential imaginary
>> (negative) values. However, upon running:
>> g_nmeig_d -f dogdic_nma.mtx -s dogdic_nma.tpr -last 3000
>> I see the following:
>> 
>> 
>> 
>> Full matrix storage format, nrow=279, ncols=279
>> Diagonalizing to find vectors 1 through 279...
>> 
>> One of the lowest 6 eigenvalues has a non-zero value.
>> 
>> This could mean that the reference structure was not
>> properly energy minimized.
>> Writing eigenvalues...
>> 
>> Back Off! I just backed up eigenval.xvg to ./#eigenval.xvg.2#
>> Writing eigenfrequencies - negative eigenvalues will be set to zero.
>> 
>> 
>> 
>> 
>> My lowest 6 eigenvalues are:
>> 1   0.0780823
>> 20.153357
>> 30.174772
>> 40.327362
>> 50.466203
>> 6 0.64693
>> 
>> 
>> I’m a little confused by the statement "One of the lowest 6 eigenvalues
>> has a non-zero value. This could mean that the reference structure was
>>not
>> properly energy minimized.” Really? Surely it out to be if one or more
>>of
>> the values is negative.  Also, there were no eigenvalues set to zero
>> (hence, I can only assume I have no negative eigenvalues).
>> 
>> Would appreciate a little insight.
>> Many thanks
>> Anthony
>> 
>> Dr Anthony Nash
>> Department of Chemistry
>> University College London
>> 
>> 
>> 
>> 
>> 
>> On 29/02/2016 16:25, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se
>>on
>> behalf of Nash, Anthony"
>> <gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
>> a.n...@ucl.ac.uk> wrote:
>> 
>>> Dear Mark and Justin,
>>> 
>>> By removing the restraints (your suggestions) it appears to have
>>>worked!
>>> 
>>> Maximum force: 9.21098e-04
>>> Writing Hessian...
>>> 
>>> This now matches the same force as my energy minimisation. Many thanks
>>>for
>>> you help.
>>> 
>>> Thanks
>>> Anthony
>>> 
>>> 
>>> 
>>> Dr Anthony Nash
>>> Department of Chemistry
>>> University College London
>>> 
>>> 
>>> 
>>> 
>>> 
>>> On 29/02/2016 14:46,
>>>"gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
>>> behalf of Mark Abraham"
>>><gromacs.org_gmx-users-boun...@maillist.sys.kth.se
>>> on behalf of mark.j.abra...@gmail.com> wrote:
>>> 
>>>> Hi,
>>>> 
>>>> An earlier mdp file suggested you were using position restraints.
>>>>There
>>>> should be no need for this, nor any problem, but what happens without
>>>> them?
>>>> 
>>>> Mark
>>>> 
>>>>> On Mon, 29 Feb 2016 15:39 Nash, Anthony <a.n...@ucl.ac.uk> wrote:
>>>>> 
>>>>> Hi Justin,
>>>>> 
>>>>> After some digging I had found that link and made some adjustments
>>>>>(as
>

Re: [gmx-users] Minimising forces for vibrational normal mode analysis

[Nash, Anthony]

Hi all,

Hoping this will be the final stumbling block. I’ve successfully minimised
my structure to a a very small max force:
Potential Energy  = -1.93736854460113e+03
Maximum force =  7.39399083846082e-04 on atom 34
Norm of force =  1.64068483698544e-04


When I perform integrator=nm there are no warning of potential imaginary
(negative) values. However, upon running:
g_nmeig_d -f dogdic_nma.mtx -s dogdic_nma.tpr -last 3000
I see the following:



Full matrix storage format, nrow=279, ncols=279
Diagonalizing to find vectors 1 through 279...

One of the lowest 6 eigenvalues has a non-zero value.

This could mean that the reference structure was not
properly energy minimized.
Writing eigenvalues...

Back Off! I just backed up eigenval.xvg to ./#eigenval.xvg.2#
Writing eigenfrequencies - negative eigenvalues will be set to zero.




My lowest 6 eigenvalues are:
1   0.0780823
 20.153357
 30.174772
 40.327362
 50.466203
 6 0.64693


I’m a little confused by the statement "One of the lowest 6 eigenvalues
has a non-zero value. This could mean that the reference structure was not
properly energy minimized.” Really? Surely it out to be if one or more of
the values is negative.  Also, there were no eigenvalues set to zero
(hence, I can only assume I have no negative eigenvalues).

Would appreciate a little insight.
Many thanks
Anthony

Dr Anthony Nash
Department of Chemistry
University College London





On 29/02/2016 16:25, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Nash, Anthony"
<gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
a.n...@ucl.ac.uk> wrote:

>Dear Mark and Justin,
>
>By removing the restraints (your suggestions) it appears to have worked!
>
>Maximum force: 9.21098e-04
>Writing Hessian...
>
>This now matches the same force as my energy minimisation. Many thanks for
>you help.
>
>Thanks
>Anthony
>
>
>
>Dr Anthony Nash
>Department of Chemistry
>University College London
>
>
>
>
>
>On 29/02/2016 14:46, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
>behalf of Mark Abraham" <gromacs.org_gmx-users-boun...@maillist.sys.kth.se
>on behalf of mark.j.abra...@gmail.com> wrote:
>
>>Hi,
>>
>>An earlier mdp file suggested you were using position restraints. There
>>should be no need for this, nor any problem, but what happens without
>>them?
>>
>>Mark
>>
>>On Mon, 29 Feb 2016 15:39 Nash, Anthony <a.n...@ucl.ac.uk> wrote:
>>
>>> Hi Justin,
>>>
>>> After some digging I had found that link and made some adjustments (as
>>> presented in the later email).
>>>
>>> After a series of energy minimisations (including switching LINCS off,
>>>and
>>> dropping the energy step to a very small number), and with the final
>>> command:
>>> grompp_d -f cg.mdp -c modic_cg_3.gro -p system.top -o modic_cg_4 -t
>>> modic_cg_3.trr
>>>
>>>
>>> Polak-Ribiere Conjugate Gradients converged to Fmax < 0.001 in 6883
>>>steps
>>> Potential Energy  = -1.92681826422996e+03
>>> Maximum force =  8.54361662283108e-04 on atom 44
>>> Norm of force =  2.71110116926712e-04
>>>
>>>
>>> However, when switching to integrator=nm, running grompp_d -f nma.mdp
>>>-c
>>> modic_cg_4.gro -p system.top -o modic_nma -t modic_cg_4.trr, and then
>>> mdrun_d, I get:
>>>
>>> Maximum force: 5.19179e+01
>>> The force is probably not small enough to ensure that you are at a
>>>minimum.
>>> Be aware that negative eigenvalues may occur
>>> when the resulting matrix is diagonalized.
>>>
>>>
>>>
>>> I¹m still struggling to yield a maximum force during the integrator=nm
>>> step, as presented from the earlier cg step. The .mdp files are
>>>identical
>>> with the exception of the integrator line.
>>>
>>> Thanks
>>> Anthony
>>>
>>>
>>> Dr Anthony Nash
>>> Department of Chemistry
>>> University College London
>>>
>>>
>>>
>>>
>>>
>>> On 29/02/2016 12:49, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se
>>>on
>>> behalf of Justin Lemkul"
>>> <gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
>>> jalem...@vt.edu> wrote:
>>>
>>> >
>>> >
>>> >On 2/29/16 3:41 AM, Nash, Anthony wrote:
>>> >> Hi Tsjerk,
>>> >>
>>

Re: [gmx-users] Minimising forces for vibrational normal mode analysis

[Nash, Anthony]

Dear Mark and Justin,

By removing the restraints (your suggestions) it appears to have worked!

Maximum force: 9.21098e-04
Writing Hessian...

This now matches the same force as my energy minimisation. Many thanks for
you help.

Thanks
Anthony



Dr Anthony Nash
Department of Chemistry
University College London





On 29/02/2016 14:46, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Mark Abraham" <gromacs.org_gmx-users-boun...@maillist.sys.kth.se
on behalf of mark.j.abra...@gmail.com> wrote:

>Hi,
>
>An earlier mdp file suggested you were using position restraints. There
>should be no need for this, nor any problem, but what happens without
>them?
>
>Mark
>
>On Mon, 29 Feb 2016 15:39 Nash, Anthony <a.n...@ucl.ac.uk> wrote:
>
>> Hi Justin,
>>
>> After some digging I had found that link and made some adjustments (as
>> presented in the later email).
>>
>> After a series of energy minimisations (including switching LINCS off,
>>and
>> dropping the energy step to a very small number), and with the final
>> command:
>> grompp_d -f cg.mdp -c modic_cg_3.gro -p system.top -o modic_cg_4 -t
>> modic_cg_3.trr
>>
>>
>> Polak-Ribiere Conjugate Gradients converged to Fmax < 0.001 in 6883
>>steps
>> Potential Energy  = -1.92681826422996e+03
>> Maximum force =  8.54361662283108e-04 on atom 44
>> Norm of force =  2.71110116926712e-04
>>
>>
>> However, when switching to integrator=nm, running grompp_d -f nma.mdp -c
>> modic_cg_4.gro -p system.top -o modic_nma -t modic_cg_4.trr, and then
>> mdrun_d, I get:
>>
>> Maximum force: 5.19179e+01
>> The force is probably not small enough to ensure that you are at a
>>minimum.
>> Be aware that negative eigenvalues may occur
>> when the resulting matrix is diagonalized.
>>
>>
>>
>> I¹m still struggling to yield a maximum force during the integrator=nm
>> step, as presented from the earlier cg step. The .mdp files are
>>identical
>> with the exception of the integrator line.
>>
>> Thanks
>> Anthony
>>
>>
>> Dr Anthony Nash
>> Department of Chemistry
>> University College London
>>
>>
>>
>>
>>
>> On 29/02/2016 12:49, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se
>>on
>> behalf of Justin Lemkul"
>> <gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
>> jalem...@vt.edu> wrote:
>>
>> >
>> >
>> >On 2/29/16 3:41 AM, Nash, Anthony wrote:
>> >> Hi Tsjerk,
>> >>
>> >> The two .mdp files are virtually identical (the only exception being
>> >>what
>> >> defines one as a conjugate-gradient, and the other for normal mode
>> >> analysis):
>> >>
>> >> CONJUGATE GRADIENT:
>> >> define = -DPOSRES
>> >> integrator  = cg
>> >> emtol   = 0.001
>> >> emstep  = 0.0002
>> >> nsteps  = 100
>> >> nstcgsteep  = 100
>> >> cutoff-scheme = verlet
>> >> nstlist = 10
>> >> ns_type = grid
>> >> rlist   = 1.4
>> >> coulombtype = PME
>> >> rcoulomb= 1.4
>> >> rvdw= 1.4
>> >> pbc = xyz
>> >>
>> >>
>> >> NORMAL MODE ANALYSIS
>> >> define = -DPOSRES
>> >> integrator = nm
>> >> emtol = 0.001
>> >> emstep = 0.0002
>> >> nsteps = 100
>> >> cutoff-scheme = verlet
>> >> nstlist = 10
>> >> ns_type = grid
>> >> rlist = 1.4
>> >> coulombtype = PME
>> >> rcoulomb = 1.4
>> >> rvdw = 1.4
>> >> pbc = xyz
>> >>
>> >>
>> >> The cg energy minimisation did NOT result in any warning about force
>>not
>> >> converging. The result I got (I just re did it now) was:
>> >> Polak-Ribiere Conjugate Gradients converged to Fmax < 0.001 in 6994
>> >>steps
>> >> Potential Energy  = -1.73087278108256e+03
>> >> Maximum force =  9.97199385029354e-04 on atom 72
>> >> Norm of force =  4.63373534634654e-04
>> >>
>> >>
>> >> But then when I run normal mode analysis (integrator=nm) I get:
>> >> Maximum force: 6.97334e+02
>> >> The force is probably not small enough to ensure that you are at a
>> >>minimum.
>> >> Be aware that negative eig

Re: [gmx-users] Minimising forces for vibrational normal mode analysis

[Nash, Anthony]

Hi Justin,

After some digging I had found that link and made some adjustments (as
presented in the later email).

After a series of energy minimisations (including switching LINCS off, and
dropping the energy step to a very small number), and with the final
command:
grompp_d -f cg.mdp -c modic_cg_3.gro -p system.top -o modic_cg_4 -t
modic_cg_3.trr


Polak-Ribiere Conjugate Gradients converged to Fmax < 0.001 in 6883 steps
Potential Energy  = -1.92681826422996e+03
Maximum force =  8.54361662283108e-04 on atom 44
Norm of force =  2.71110116926712e-04


However, when switching to integrator=nm, running grompp_d -f nma.mdp -c
modic_cg_4.gro -p system.top -o modic_nma -t modic_cg_4.trr, and then
mdrun_d, I get:

Maximum force: 5.19179e+01
The force is probably not small enough to ensure that you are at a minimum.
Be aware that negative eigenvalues may occur
when the resulting matrix is diagonalized.



I¹m still struggling to yield a maximum force during the integrator=nm
step, as presented from the earlier cg step. The .mdp files are identical
with the exception of the integrator line.

Thanks
Anthony 


Dr Anthony Nash
Department of Chemistry
University College London





On 29/02/2016 12:49, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Justin Lemkul"
<gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
jalem...@vt.edu> wrote:

>
>
>On 2/29/16 3:41 AM, Nash, Anthony wrote:
>> Hi Tsjerk,
>>
>> The two .mdp files are virtually identical (the only exception being
>>what
>> defines one as a conjugate-gradient, and the other for normal mode
>> analysis):
>>
>> CONJUGATE GRADIENT:
>> define = -DPOSRES
>> integrator  = cg
>> emtol   = 0.001
>> emstep  = 0.0002
>> nsteps  = 100
>> nstcgsteep  = 100
>> cutoff-scheme = verlet
>> nstlist = 10
>> ns_type = grid
>> rlist   = 1.4
>> coulombtype = PME
>> rcoulomb= 1.4
>> rvdw= 1.4
>> pbc = xyz
>>
>>
>> NORMAL MODE ANALYSIS
>> define = -DPOSRES
>> integrator = nm
>> emtol = 0.001
>> emstep = 0.0002
>> nsteps = 100
>> cutoff-scheme = verlet
>> nstlist = 10
>> ns_type = grid
>> rlist = 1.4
>> coulombtype = PME
>> rcoulomb = 1.4
>> rvdw = 1.4
>> pbc = xyz
>>
>>
>> The cg energy minimisation did NOT result in any warning about force not
>> converging. The result I got (I just re did it now) was:
>> Polak-Ribiere Conjugate Gradients converged to Fmax < 0.001 in 6994
>>steps
>> Potential Energy  = -1.73087278108256e+03
>> Maximum force =  9.97199385029354e-04 on atom 72
>> Norm of force =  4.63373534634654e-04
>>
>>
>> But then when I run normal mode analysis (integrator=nm) I get:
>> Maximum force: 6.97334e+02
>> The force is probably not small enough to ensure that you are at a
>>minimum.
>> Be aware that negative eigenvalues may occur
>> when the resulting matrix is diagonalized.
>>
>>
>> My work flow:
>> grompp_d -f cg.mdp -c modic_en.gro -p system.top -o modic_cg
>> mdrun_d -deffnm modic_cg
>> grompp_d -f nma.mdp -c modic_cg.gro -p system.top -o modic_nma
>> mdrun_d -deffnm modic_nma
>>
>
>Relying only on .gro format is insufficient precision to continue from
>minimization to NMA.  Make sure to pass the .trr to grompp -t.
>
>http://www.gromacs.org/Documentation/How-tos/Normal_Mode_Analysis
>
>-Justin
>
>-- 
>==
>
>Justin A. Lemkul, Ph.D.
>Ruth L. Kirschstein NRSA Postdoctoral Fellow
>
>Department of Pharmaceutical Sciences
>School of Pharmacy
>Health Sciences Facility II, Room 629
>University of Maryland, Baltimore
>20 Penn St.
>Baltimore, MD 21201
>
>jalem...@outerbanks.umaryland.edu | (410) 706-7441
>http://mackerell.umaryland.edu/~jalemkul
>
>==
>-- 
>Gromacs Users mailing list
>
>* Please search the archive at
>http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>posting!
>
>* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
>* For (un)subscribe requests visit
>https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>send a mail to gmx-users-requ...@gromacs.org.

-- 
Gromacs Users mailing list

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

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

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

Re: [gmx-users] Minimising forces for vibrational normal mode analysis

[Nash, Anthony]

I noticed that I hadn’t included LINCS parameters in the mdp file. I have
included and currently running using:

continuation = no
constraints = all-bonds ; all bonds (even heavy atom-H
bonds) constrained
lincs_iter  = 2 ; accuracy of LINCS
lincs_order = 8 ; also related to accuracy


I upped lincs_iter from 1 to 2, and lincs_order from 4 to 8. The manual
says this is ideal for increasing accuracy during energy minimisation.

After performing this run I’m now getting forces beyond Fmax and I’ve
reached the maximum number of steps.
Potential Energy  = -1.71320900832000e+03
Maximum force =  1.80888456004520e+01 on atom 13
Norm of force =  7.70125465905422e+00


Prior to including LINCS, the system was finishing within 1000 steps, so I
feel as though I’m heading down the right paths in trying to get these
forces down first now that I’ve included lincs.




Dr Anthony Nash
Department of Chemistry
University College London





On 29/02/2016 08:41, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Nash, Anthony"
<gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
a.n...@ucl.ac.uk> wrote:

>Hi Tsjerk,
>
>The two .mdp files are virtually identical (the only exception being what
>defines one as a conjugate-gradient, and the other for normal mode
>analysis):
>
>CONJUGATE GRADIENT:
>define = -DPOSRES
>integrator  = cg
>emtol   = 0.001
>emstep  = 0.0002
>nsteps  = 100
>nstcgsteep  = 100
>cutoff-scheme = verlet
>nstlist = 10
>ns_type = grid
>rlist   = 1.4
>coulombtype = PME
>rcoulomb= 1.4
>rvdw= 1.4
>pbc = xyz
>
>
>NORMAL MODE ANALYSIS
>define = -DPOSRES
>integrator = nm 
>emtol = 0.001 
>emstep = 0.0002 
>nsteps = 100 
>cutoff-scheme = verlet
>nstlist = 10
>ns_type = grid 
>rlist = 1.4 
>coulombtype = PME 
>rcoulomb = 1.4 
>rvdw = 1.4 
>pbc = xyz 
>
>
>The cg energy minimisation did NOT result in any warning about force not
>converging. The result I got (I just re did it now) was:
>Polak-Ribiere Conjugate Gradients converged to Fmax < 0.001 in 6994 steps
>Potential Energy  = -1.73087278108256e+03
>Maximum force =  9.97199385029354e-04 on atom 72
>Norm of force =  4.63373534634654e-04
>
>
>But then when I run normal mode analysis (integrator=nm) I get:
>Maximum force: 6.97334e+02
>The force is probably not small enough to ensure that you are at a
>minimum.
>Be aware that negative eigenvalues may occur
>when the resulting matrix is diagonalized.
>
>
>My work flow:
>grompp_d -f cg.mdp -c modic_en.gro -p system.top -o modic_cg
>mdrun_d -deffnm modic_cg
>grompp_d -f nma.mdp -c modic_cg.gro -p system.top -o modic_nma
>mdrun_d -deffnm modic_nma
>
>
>
>
>
>
>I sincerely hope I’m doing something wrong. That way it’ll be easier to
>solve. 
>
>Many thanks
>Anthony
>
>
>Dr Anthony Nash
>Department of Chemistry
>University College London
>
>
>
>
>
>On 29/02/2016 08:03, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
>behalf of Tsjerk Wassenaar"
><gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
>tsje...@gmail.com> wrote:
>
>>Hi Anthony,
>>
>>You can paste your exact workflow commands and mdp files, but that would
>>be
>>to just check whether you overlooked something. If you assert that the
>>mdp
>>files are the same, except for the integrator and you really did use the
>>final minimized structure as input for the nm run, then it seems
>>something's fishy. You could make a run input file with integrator=cg and
>>integrator=nm and compare the two tpr files to see if something was
>>changed
>>implicitly.
>>
>>Cheers,
>>
>>Tsjerk
>>
>>On Mon, Feb 29, 2016 at 6:57 AM, Nash, Anthony <a.n...@ucl.ac.uk> wrote:
>>
>>> Hi Tsjerk,
>>>
>>> Compiled in double precision. I’ve alternated between steepest and
>>> conjugate gradient. I am slightly confused over why the energy
>>> minimisation reports a "Maximum force =  3.72351263315387e-05 on
>>>atom
>>> 34” only then for integrator=nm to report "Maximum force: 6.02183e+02”,
>>> which will presumably return in an imaginary frequency (I haven’t tried
>>>-
>>> yet even if it didn’t I would be suspicious given the earlier warning).
>>>
>>> Thanks
>>> Anthony
>>>
>>>
>>> Dr Anthony Nash
>>> Department of Chemistry
>>> University College London
>>>
>>>
>>&g

Re: [gmx-users] Minimising forces for vibrational normal mode analysis

[Nash, Anthony]

Hi Tsjerk,

The two .mdp files are virtually identical (the only exception being what
defines one as a conjugate-gradient, and the other for normal mode
analysis):

CONJUGATE GRADIENT:
define = -DPOSRES
integrator  = cg
emtol   = 0.001
emstep  = 0.0002
nsteps  = 100
nstcgsteep  = 100
cutoff-scheme = verlet
nstlist = 10
ns_type = grid
rlist   = 1.4
coulombtype = PME
rcoulomb= 1.4
rvdw= 1.4
pbc = xyz


NORMAL MODE ANALYSIS
define = -DPOSRES
integrator = nm 
emtol = 0.001 
emstep = 0.0002 
nsteps = 100 
cutoff-scheme = verlet
nstlist = 10
ns_type = grid 
rlist = 1.4 
coulombtype = PME 
rcoulomb = 1.4 
rvdw = 1.4 
pbc = xyz 


The cg energy minimisation did NOT result in any warning about force not
converging. The result I got (I just re did it now) was:
Polak-Ribiere Conjugate Gradients converged to Fmax < 0.001 in 6994 steps
Potential Energy  = -1.73087278108256e+03
Maximum force =  9.97199385029354e-04 on atom 72
Norm of force =  4.63373534634654e-04


But then when I run normal mode analysis (integrator=nm) I get:
Maximum force: 6.97334e+02
The force is probably not small enough to ensure that you are at a minimum.
Be aware that negative eigenvalues may occur
when the resulting matrix is diagonalized.


My work flow:
grompp_d -f cg.mdp -c modic_en.gro -p system.top -o modic_cg
mdrun_d -deffnm modic_cg
grompp_d -f nma.mdp -c modic_cg.gro -p system.top -o modic_nma
mdrun_d -deffnm modic_nma






I sincerely hope I’m doing something wrong. That way it’ll be easier to
solve. 

Many thanks
Anthony


Dr Anthony Nash
Department of Chemistry
University College London





On 29/02/2016 08:03, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Tsjerk Wassenaar"
<gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
tsje...@gmail.com> wrote:

>Hi Anthony,
>
>You can paste your exact workflow commands and mdp files, but that would
>be
>to just check whether you overlooked something. If you assert that the mdp
>files are the same, except for the integrator and you really did use the
>final minimized structure as input for the nm run, then it seems
>something's fishy. You could make a run input file with integrator=cg and
>integrator=nm and compare the two tpr files to see if something was
>changed
>implicitly.
>
>Cheers,
>
>Tsjerk
>
>On Mon, Feb 29, 2016 at 6:57 AM, Nash, Anthony <a.n...@ucl.ac.uk> wrote:
>
>> Hi Tsjerk,
>>
>> Compiled in double precision. I’ve alternated between steepest and
>> conjugate gradient. I am slightly confused over why the energy
>> minimisation reports a "Maximum force =  3.72351263315387e-05 on
>>atom
>> 34” only then for integrator=nm to report "Maximum force: 6.02183e+02”,
>> which will presumably return in an imaginary frequency (I haven’t tried
>>-
>> yet even if it didn’t I would be suspicious given the earlier warning).
>>
>> Thanks
>> Anthony
>>
>>
>> Dr Anthony Nash
>> Department of Chemistry
>> University College London
>>
>>
>>
>>
>>
>> On 28/02/2016 22:04, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se
>>on
>> behalf of Tsjerk Wassenaar"
>> <gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
>> tsje...@gmail.com> wrote:
>>
>> >Hi Anthony,
>> >
>> >You do not state whether Gromacs was compiled in double precision or
>>not.
>> >It should be for this stuff. In addition, you can try another
>>minimization
>> >method. Sometimes alternating minimization methods may help to reach a
>> >proper minimum.
>> >
>> >Cheers,
>> >
>> >Tsjerk
>> >
>> >On Sun, Feb 28, 2016 at 11:27 AM, Nash, Anthony <a.n...@ucl.ac.uk>
>>wrote:
>> >
>> >> Hi all,
>> >>
>> >> I would like to pull out the vibrational normal modes using gromacs
>> >>over a
>> >> customised fragment to compare back with the original QM frequency
>> >> analysis.
>> >>
>> >> I¹ve performed an integrator=cg over my structure, and monitored the
>> >> potential energy which converges. The forces also converge beneath
>>the
>> >> requested precision (as 0.0001, as per gromacs manual). The message:
>> >>
>> >> 
>> >>
>> >> Polak-Ribiere Conjugate Gradients:
>> >>Tolerance (Fmax)   =  1.0e-04
>> >>Number of steps=  100
>> >>F-max =  7.35939e+03 on atom 79
>> >>F-Norm=  1.68505e+03

Re: [gmx-users] Minimising forces for vibrational normal mode analysis

[Nash, Anthony]

Hi Tsjerk,

Compiled in double precision. I’ve alternated between steepest and
conjugate gradient. I am slightly confused over why the energy
minimisation reports a "Maximum force =  3.72351263315387e-05 on atom
34” only then for integrator=nm to report "Maximum force: 6.02183e+02”,
which will presumably return in an imaginary frequency (I haven’t tried -
yet even if it didn’t I would be suspicious given the earlier warning).

Thanks
Anthony 


Dr Anthony Nash
Department of Chemistry
University College London





On 28/02/2016 22:04, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Tsjerk Wassenaar"
<gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
tsje...@gmail.com> wrote:

>Hi Anthony,
>
>You do not state whether Gromacs was compiled in double precision or not.
>It should be for this stuff. In addition, you can try another minimization
>method. Sometimes alternating minimization methods may help to reach a
>proper minimum.
>
>Cheers,
>
>Tsjerk
>
>On Sun, Feb 28, 2016 at 11:27 AM, Nash, Anthony <a.n...@ucl.ac.uk> wrote:
>
>> Hi all,
>>
>> I would like to pull out the vibrational normal modes using gromacs
>>over a
>> customised fragment to compare back with the original QM frequency
>> analysis.
>>
>> I¹ve performed an integrator=cg over my structure, and monitored the
>> potential energy which converges. The forces also converge beneath the
>> requested precision (as 0.0001, as per gromacs manual). The message:
>>
>> 
>>
>> Polak-Ribiere Conjugate Gradients:
>>Tolerance (Fmax)   =  1.0e-04
>>Number of steps=  100
>>F-max =  7.35939e+03 on atom 79
>>F-Norm=  1.68505e+03
>>
>> writing lowest energy coordinates.
>>
>> Polak-Ribiere Conjugate Gradients converged to Fmax < 0.0001 in 8350
>>steps
>> Potential Energy  = -1.67001242395943e+03
>> Maximum force =  3.72351263315387e-05 on atom 34
>> Norm of force =  1.20696019277311e-05
>>
>>
>> 
>>
>>
>> When I then come to run integrator=nm I get a maximum force at odds with
>> the maximum force reported at the final stage of my energy minimisation:
>>
>> 
>>
>> Maximum force: 6.02183e+02
>> The force is probably not small enough to ensure that you are at a
>>minimum.
>> Be aware that negative eigenvalues may occur
>> when the resulting matrix is diagonalized.
>>
>> 
>>
>>
>>
>> When I decrease the force tolerance in the integrator=cg energy
>> minimisation to 0.1 I end up with a poorer force convergence
>>(although
>> the potential energy is almost the same, but also the integrator=nm will
>> result in the same measure of maximum force):
>>
>> 
>>
>> Polak-Ribiere Conjugate Gradients:
>>Tolerance (Fmax)   =  1.0e-05
>>Number of steps=  100
>>F-max =  7.35939e+03 on atom 79
>>F-Norm=  1.68505e+03
>>
>> Energy minimization has stopped, but the forces have not converged to
>>the
>> requested precision Fmax < 1e-05 (which may not be possible for your
>> system).
>> It stopped because the algorithm tried to make a new step whose size was
>> too
>> small, or there was no change in the energy since last step. Either
>>way, we
>> regard the minimization as converged to within the available machine
>> precision, given your starting configuration and EM parameters.
>>
>> writing lowest energy coordinates.
>>
>> Polak-Ribiere Conjugate Gradients converged to machine precision in 9839
>> steps,
>> but did not reach the requested Fmax < 1e-05.
>> Potential Energy  = -1.67001242392665e+03
>> Maximum force =  1.63690887039393e-03 on atom 34
>> Norm of force =  3.34598700357485e-04
>> 
>>
>>
>>
>> I don¹t want to end up with any imaginary values when I diagonalise the
>> hessian, any idea how to improve this performance? I am concerned with
>>the
>> output "Energy minimization has stopped, but the forces have not
>>converged
>> to the requested precision Fmax < 1e-05 (which may not be possible for
>> your system).² a a possible indication to the computational limitation
>>of
>> my machine.
>>
>> Many thanks
>> Anthony
>>
>>
>> --
>> Gromacs Users mailing list
>>
>> * Please search the archive

[gmx-users] Minimising forces for vibrational normal mode analysis

[Nash, Anthony]

Hi all,

I would like to pull out the vibrational normal modes using gromacs over a
customised fragment to compare back with the original QM frequency
analysis. 

I¹ve performed an integrator=cg over my structure, and monitored the
potential energy which converges. The forces also converge beneath the
requested precision (as 0.0001, as per gromacs manual). The message:



Polak-Ribiere Conjugate Gradients:
   Tolerance (Fmax)   =  1.0e-04
   Number of steps=  100
   F-max =  7.35939e+03 on atom 79
   F-Norm=  1.68505e+03

writing lowest energy coordinates.

Polak-Ribiere Conjugate Gradients converged to Fmax < 0.0001 in 8350 steps
Potential Energy  = -1.67001242395943e+03
Maximum force =  3.72351263315387e-05 on atom 34
Norm of force =  1.20696019277311e-05





When I then come to run integrator=nm I get a maximum force at odds with
the maximum force reported at the final stage of my energy minimisation:



Maximum force: 6.02183e+02
The force is probably not small enough to ensure that you are at a minimum.
Be aware that negative eigenvalues may occur
when the resulting matrix is diagonalized.




 
When I decrease the force tolerance in the integrator=cg energy
minimisation to 0.1 I end up with a poorer force convergence (although
the potential energy is almost the same, but also the integrator=nm will
result in the same measure of maximum force):



Polak-Ribiere Conjugate Gradients:
   Tolerance (Fmax)   =  1.0e-05
   Number of steps=  100
   F-max =  7.35939e+03 on atom 79
   F-Norm=  1.68505e+03

Energy minimization has stopped, but the forces have not converged to the
requested precision Fmax < 1e-05 (which may not be possible for your
system).
It stopped because the algorithm tried to make a new step whose size was
too
small, or there was no change in the energy since last step. Either way, we
regard the minimization as converged to within the available machine
precision, given your starting configuration and EM parameters.

writing lowest energy coordinates.

Polak-Ribiere Conjugate Gradients converged to machine precision in 9839
steps,
but did not reach the requested Fmax < 1e-05.
Potential Energy  = -1.67001242392665e+03
Maximum force =  1.63690887039393e-03 on atom 34
Norm of force =  3.34598700357485e-04




I don¹t want to end up with any imaginary values when I diagonalise the
hessian, any idea how to improve this performance? I am concerned with the
output "Energy minimization has stopped, but the forces have not converged
to the requested precision Fmax < 1e-05 (which may not be possible for
your system).² a a possible indication to the computational limitation of
my machine. 

Many thanks
Anthony 


-- 
Gromacs Users mailing list

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

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

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

Re: [gmx-users] generating an initial structure with a gromacs compatible tool

[Nash, Anthony]

Hi Mark,

I’ve double checked and triple checked the output from Avogadro. It does
appear (within the constraints of my version and installation) that
resaving an already correct .pdb format will shift the X coordinate left
by one space. I’ve reported this issue.

Thank you for that information regarding the “renumbering”, this should
save me a lot of time.

With regards to my original problem, I fixed it this morning. I took your
recommendation to triple the atom ordering in my new residue and that’s
when I noticed the awful mistake! Rather than using the atom names in the
[ bonds ] entry of aminoacids.rtp, I had used their atom index within the
residue. Suffice to say grompp continued to work, but when I corrected
this mistake I immediately noticed an additional 50 bonds had been
processed by grompp. I’ve performed em, NVT and NPT. My structure is very
stable :) 

Many thanks for helping me trouble shoot.

Anthony




On 25/02/2016 13:12, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Mark Abraham" <gromacs.org_gmx-users-boun...@maillist.sys.kth.se
on behalf of mark.j.abra...@gmail.com> wrote:

>Hi,
>
>On Wed, Feb 24, 2016 at 4:00 PM Nash, Anthony <a.n...@ucl.ac.uk> wrote:
>
>> Hi Mark,
>>
>> I’m afraid I am not sure what you mean by your final point "Thus,
>>avoiding
>> the issue by not doing reordering.”
>>
>> I don’t think I’ve been as clear as I could have been, sorry about
>>that. I
>> am only using avogadro to attach the two glycine residues, NGLY and CGLY
>> to the backbone of my brand new fragment whose geometry comes from the
>> output of a HF/6-31G(p) optimisation. I load the Gaussian .pdb into
>> avogadro, attach the two glycine residues, and then resave without
>> disturbing the coordinates of the original fragment (of which the eq
>>bond
>> length and angles are based on). Unfortunately, going back to my
>>original
>> question of “does gromacs have its own editing tool”, avogadro adjusts
>>the
>> formats of the original Gaussian .pdb file (including deviating the
>> X-coordinate entry by a single character shift to the left, causing
>> pdb2gmx to throw out a warning for every atom in the system)
>
>
>That sounds like a bug in Avogadro, and you should complain about it. pdb
>is a fixed-column format.
>
>
>> plus puts in
>> the terminal glycines at the end of the file. This is no good as
>>pdb2gmx’s
>> understanding of amber requires the termini of each chain to be defined
>>by
>> the presence of NGLY-FRAGMENT-GCLY.
>>
>
>But this has a trivial fix just copy the NGLY fragment and put it at
>the start. The numbering doesn't matter, so long as the numbers do change
>"reasonably."
>
>Either way, I’ve successfully walked through the complete steps of
>> Gaussian optimisation and RESP derivation to a running NPT gromacs
>> simulation for an almost identical peptide residue. The fact that the
>> second one isn’t working is probably indicative of my atom ordering
>>being
>> skewered - perhaps by just two atoms!
>>
>
>Well, make a vacuum system with just ngly-fragment-cgly and do gmx dump on
>the .tpr. There's no docs on how that format works, but you can see what
>grompp has made of your topology, which can function as that "second pair
>of eyes." Just be alert that all the indexing in the dump starts from 0
>(because that's convenient in C programs).
>
>Mark
>
>Thanks
>> Anthony
>>
>> Dr Anthony Nash
>> Department of Chemistry
>> University College London
>>
>>
>>
>>
>>
>> On 24/02/2016 14:33, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se
>>on
>> behalf of Mark Abraham"
>><gromacs.org_gmx-users-boun...@maillist.sys.kth.se
>> on behalf of mark.j.abra...@gmail.com> wrote:
>>
>> >Hi,
>> >
>> >On Wed, Feb 24, 2016 at 3:26 PM Nash, Anthony <a.n...@ucl.ac.uk> wrote:
>> >
>> >> Hi Mark,
>> >>
>> >> When you generate a peptide sequences in Avogadro the atom name
>>order in
>> >> the .pdb for NGLY (which is just GLY and required renaming) is
>> >>NHCCOH
>> >> (if my memory serves me right - ignore the shortening of the names),
>> >> whilst in the .rtp file of amberffsb.99 it is NHHHCHHCO. I’ve always
>> >> reordered by hand. You mention that atom order isn’t significant
>>until
>> >> grompp, are you suggesting the pdb2gmx will understand *any* order of
>> >> atoms within a residue?
>> >>
>> >
>> >Try it :-) You have a working case, so swap the o

Re: [gmx-users] generating an initial structure with a gromacs compatible tool

[Nash, Anthony]

Hi Mark,

I’m afraid I am not sure what you mean by your final point "Thus, avoiding
the issue by not doing reordering.”

I don’t think I’ve been as clear as I could have been, sorry about that. I
am only using avogadro to attach the two glycine residues, NGLY and CGLY
to the backbone of my brand new fragment whose geometry comes from the
output of a HF/6-31G(p) optimisation. I load the Gaussian .pdb into
avogadro, attach the two glycine residues, and then resave without
disturbing the coordinates of the original fragment (of which the eq bond
length and angles are based on). Unfortunately, going back to my original
question of “does gromacs have its own editing tool”, avogadro adjusts the
formats of the original Gaussian .pdb file (including deviating the
X-coordinate entry by a single character shift to the left, causing
pdb2gmx to throw out a warning for every atom in the system) plus puts in
the terminal glycines at the end of the file. This is no good as pdb2gmx’s
understanding of amber requires the termini of each chain to be defined by
the presence of NGLY-FRAGMENT-GCLY.

Either way, I’ve successfully walked through the complete steps of
Gaussian optimisation and RESP derivation to a running NPT gromacs
simulation for an almost identical peptide residue. The fact that the
second one isn’t working is probably indicative of my atom ordering being
skewered - perhaps by just two atoms!

Thanks
Anthony 

Dr Anthony Nash
Department of Chemistry
University College London





On 24/02/2016 14:33, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Mark Abraham" <gromacs.org_gmx-users-boun...@maillist.sys.kth.se
on behalf of mark.j.abra...@gmail.com> wrote:

>Hi,
>
>On Wed, Feb 24, 2016 at 3:26 PM Nash, Anthony <a.n...@ucl.ac.uk> wrote:
>
>> Hi Mark,
>>
>> When you generate a peptide sequences in Avogadro the atom name order in
>> the .pdb for NGLY (which is just GLY and required renaming) is
>>NHCCOH
>> (if my memory serves me right - ignore the shortening of the names),
>> whilst in the .rtp file of amberffsb.99 it is NHHHCHHCO. I’ve always
>> reordered by hand. You mention that atom order isn’t significant until
>> grompp, are you suggesting the pdb2gmx will understand *any* order of
>> atoms within a residue?
>>
>
>Try it :-) You have a working case, so swap the order of the atoms in the
>pdb2gmx input coordinate file, and re-run your scripts.
>
>With regards to your point about constraints - yes, you are right. I’m
>> using an energy minimisation .mdp input file without any mention of
>> constraints. I’ve never once had to do this within an energy
>>minimisation
>> (I must have very fortunate starting structures up until yesterday)
>>until
>> I come to use NVT and/or NPT steeping. I’ve just put on LINCS and reran
>> the energy minimiation, see below for the output:
>>
>
>You aren't going to be able to use what Avogadro produces without
>understanding what it is producing. Just turning on constraints isn't
>going
>to result in a valid model if it was "parameterized" for something else.
>Clearly the output coordinate file is not very close to the constrained
>geometry, but you need to understand why that is before you know how to
>handle it.
>
>
>> Regarding the atom ordering - this is something I suspect, but I’ve
>> trawled through every gromacs file I’ve changed, and input geometry and
>>it
>> all seems to be aligned. I suspect I need a fresh set of eyes.
>>
>
>Thus, avoiding the issue by not doing reordering.
>
>Mark
>
>
>> Thanks
>> Anthony
>>
>>
>> --
>> Steepest Descents:
>>Tolerance (Fmax)   =  1.0e+01
>>Number of steps=   20
>>
>>
>> Step 40, time 0.04 (ps)  LINCS WARNING
>> relative constraint deviation after LINCS:
>> rms 0.007405, max 0.034928 (between atoms 73 and 75)
>> bonds that rotated more than 30 degrees:
>>  atom 1 atom 2  angle  previous, current, constraint length
>>  71 72   31.40.1014   0.1006  0.1010
>>
>>
>> Energy minimization has stopped, but the forces have not converged to
>>the
>> requested precision Fmax < 10 (which may not be possible for your
>>system).
>> It
>> stopped because the algorithm tried to make a new step whose size was
>>too
>> small, or there was no change in the energy since last step. Either
>>way, we
>> regard the minimization as converged to within the available machine
>> precision, given your starting configuration and EM parameters.
>> You might need to increase your constraint accuracy, or turn
>> off

Re: [gmx-users] generating an initial structure with a gromacs compatible tool

[Nash, Anthony]

Hi Mark,

When you generate a peptide sequences in Avogadro the atom name order in
the .pdb for NGLY (which is just GLY and required renaming) is NHCCOH
(if my memory serves me right - ignore the shortening of the names),
whilst in the .rtp file of amberffsb.99 it is NHHHCHHCO. I’ve always
reordered by hand. You mention that atom order isn’t significant until
grompp, are you suggesting the pdb2gmx will understand *any* order of
atoms within a residue?

With regards to your point about constraints - yes, you are right. I’m
using an energy minimisation .mdp input file without any mention of
constraints. I’ve never once had to do this within an energy minimisation
(I must have very fortunate starting structures up until yesterday) until
I come to use NVT and/or NPT steeping. I’ve just put on LINCS and reran
the energy minimiation, see below for the output:

Regarding the atom ordering - this is something I suspect, but I’ve
trawled through every gromacs file I’ve changed, and input geometry and it
all seems to be aligned. I suspect I need a fresh set of eyes.

Thanks
Anthony 


--
Steepest Descents:
   Tolerance (Fmax)   =  1.0e+01
   Number of steps=   20


Step 40, time 0.04 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 0.007405, max 0.034928 (between atoms 73 and 75)
bonds that rotated more than 30 degrees:
 atom 1 atom 2  angle  previous, current, constraint length
 71 72   31.40.1014   0.1006  0.1010


Energy minimization has stopped, but the forces have not converged to the
requested precision Fmax < 10 (which may not be possible for your system).
It
stopped because the algorithm tried to make a new step whose size was too
small, or there was no change in the energy since last step. Either way, we
regard the minimization as converged to within the available machine
precision, given your starting configuration and EM parameters.
You might need to increase your constraint accuracy, or turn
off constraints altogether (set constraints = none in mdp file)
―






Dr Anthony Nash
Department of Chemistry
University College London





On 24/02/2016 14:04, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Mark Abraham" <gromacs.org_gmx-users-boun...@maillist.sys.kth.se
on behalf of mark.j.abra...@gmail.com> wrote:

>Hi,
>
>On Tue, Feb 23, 2016 at 11:03 AM Nash, Anthony <a.n...@ucl.ac.uk> wrote:
>
>> Hi all,
>> Is there a friendly Gromacs compatible tool for generating a .gro/.pdb
>> file using a specific forcefield topology specification within Gromacs
>> itself? For context:
>>
>> I¹m in the process of fully parameterising five custom protein residues
>> for the amber forcefield from ab initio calculations. Fragment #1 has
>>been
>> very successfully, with a production NPT simulation showing very little
>> deviation from the equilibrium ab initio structure.
>>
>> Fragment #2, on the other hand, is driving me crazy. I derive partial
>> charges, force constants, then I take the equilibrium structure, and
>> replace the capped backbone ends (required for RESP) with NGLY and
>>CGLY. I
>> am using Avogadro to do this, which spits out a awful .pdb file which
>> requires a lot of rearranging to be compatible with the atom order of
>>the
>> residues inside aminoacid.rtp in the amber99sb.ff.
>
>
>But the .rtp entry is based on atom names. Atom order isn't significant
>until grompp (and then it merely warns if the atom order of the input .gro
>file does not match that of the .top). And you have an .rtp entry already,
>what are you getting from Avogadro anyway?
>
>I hold the eq_structure
>> fixed in avogadro and run an energy equilibrium over the new adjoining
>> glycine residues. All appears fine, until I run a Gromacs energy
>> minimisation. The entire customer residue expands gently (beyond the eq
>> distances, by a few angstrom),
>
>
>Clearly your energy minimization isn't using constraints. Should it? You
>need to understand the basis under which Avogadro is parameterizing for
>this force field...
>
>
>> and then an NVT simulation throws out a lot
>> of LINCS warnings.
>>
>> I am pretty confident that this is the fault of my initial structure, in
>> particular the angles (I¹m going to double check the FCs).
>
>
>The simplest explanation is that you are mangling the atom order so that
>the topology mdrun understands isn't the one you intend. But that comes
>down to what you have in your .rtp/.itp files.
>
>
>> But mean while
>> using Avogadro to build test structures is taking forever! Any
>> alternatives?
>
>
>acpype is popular conversion tool for AMBER topologies, I 

[gmx-users] structure expanding beyond eq values

[Nash, Anthony]

Hi all,

Any thoughts on what could be causing my structure to expand and distort
well beyond (about 2 to 3 angstrom with some distorted angles) the
equilibrium bond lengths during energy minimisation?

I¹ve fully parameterised two new fragments, which include new atom types
and force constants. The first fragment (sitting between two glycine
residues) behaves very well. The second fragment (separate simulation
altogether, similar set up, only difference is an added ethane group to
the central ring of the fragment) doesn¹t get past the energy minimum
stage.

As an example on bond length force constants:

Fragment #1
MC7 MN4 1 0.146 86531.610

Fragment #2
GC11 GN3 1 0.145 57828.550


This bond is from the same part of a circular structure, the only
difference is a functional group elsewhere in the ring. That particular
bond in fragment #2 goes from 0.145 to 0.3 angstrom. I¹ve tried putting
the initial structure, and the broken-minimised structure, into an NVT but
it results in LINCS errors.

Note: the initial structure came from a HF/6-31G(d) optimisation.

Many thanks
Anthony


Dr Anthony Nash
Department of Chemistry
University College London


-- 
Gromacs Users mailing list

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

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

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

[gmx-users] generating an initial structure with a gromacs compatible tool

[Nash, Anthony]

Hi all,
Is there a friendly Gromacs compatible tool for generating a .gro/.pdb
file using a specific forcefield topology specification within Gromacs
itself? For context:

I¹m in the process of fully parameterising five custom protein residues
for the amber forcefield from ab initio calculations. Fragment #1 has been
very successfully, with a production NPT simulation showing very little
deviation from the equilibrium ab initio structure.

Fragment #2, on the other hand, is driving me crazy. I derive partial
charges, force constants, then I take the equilibrium structure, and
replace the capped backbone ends (required for RESP) with NGLY and CGLY. I
am using Avogadro to do this, which spits out a awful .pdb file which
requires a lot of rearranging to be compatible with the atom order of the
residues inside aminoacid.rtp in the amber99sb.ff. I hold the eq_structure
fixed in avogadro and run an energy equilibrium over the new adjoining
glycine residues. All appears fine, until I run a Gromacs energy
minimisation. The entire customer residue expands gently (beyond the eq
distances, by a few angstrom), and then an NVT simulation throws out a lot
of LINCS warnings. 

I am pretty confident that this is the fault of my initial structure, in
particular the angles (I¹m going to double check the FCs). But mean while
using Avogadro to build test structures is taking forever! Any
alternatives?

Many thanks
Anthony

Dr Anthony Nash
Department of Chemistry
University College London

-- 
Gromacs Users mailing list

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

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

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

Re: [gmx-users] Unknown bond_atomtype

[Nash, Anthony]

Thanks Justin,

I think that¹s everything I need to know.

Kind regards
Anthony

Dr Anthony Nash
Department of Chemistry
University College London





On 17/02/2016 22:18, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Justin Lemkul"
<gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
jalem...@vt.edu> wrote:

>
>
>On 2/17/16 5:11 PM, Nash, Anthony wrote:
>> Hi,
>>
>> Thanks for the information, Mark. I fully parameterised the bonded terms
>> using QM data and some software I wrote. Now that I¹ve sorted out
>> ffnonbonded [atomtypes], all that remains are some missing angle terms.
>>I
>> want to have the forcefield files correct so an MD simulation can be
>> generated with very little changes to the topol file after running
>> pdb2gmx. I do have a couple of questions that I hope you could help me
>> with.
>>
>> 1) For grompp to compile successfully, must every molecule/residue have
>>a
>> complete set of dihedral angle force constants?
>
>grompp will complain about any missing parameters, triggering a failure.
>
>> 2) It appears as though the residue to residue bond connectivity in my
>> force field files aren¹t correctly parameterised. In NGLY-MOD-CGLY, a
>> bonded term is added to the topol for the C in NGLY to the N in CGLY. To
>> avoid this must I add an entry such as -C N in the [ bonds ] section of
>> the [ MOD1 ] residue definition in aminoacids.rtp? I¹m unsure what ³-C²
>> means.
>>
>
>- indicates an atom in the previous residue, + means in the next residue.
> So -C 
>is the C in the previous residue.  If you're getting bonds between NGLY
>and 
>CGLY, ignoring the intervening residue, you probably either haven't set
>your 
>modified residue as Protein in residuetypes.dat or it doesn't have a C
>(which it 
>should, if it's an amino acid.
>
>-Justin
>
>-- 
>==
>
>Justin A. Lemkul, Ph.D.
>Ruth L. Kirschstein NRSA Postdoctoral Fellow
>
>Department of Pharmaceutical Sciences
>School of Pharmacy
>Health Sciences Facility II, Room 629
>University of Maryland, Baltimore
>20 Penn St.
>Baltimore, MD 21201
>
>jalem...@outerbanks.umaryland.edu | (410) 706-7441
>http://mackerell.umaryland.edu/~jalemkul
>
>==
>-- 
>Gromacs Users mailing list
>
>* Please search the archive at
>http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>posting!
>
>* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
>* For (un)subscribe requests visit
>https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>send a mail to gmx-users-requ...@gromacs.org.

-- 
Gromacs Users mailing list

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

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

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

Re: [gmx-users] Unknown bond_atomtype

[Nash, Anthony]

Hi,

Thanks for the information, Mark. I fully parameterised the bonded terms
using QM data and some software I wrote. Now that I’ve sorted out
ffnonbonded [atomtypes], all that remains are some missing angle terms. I
want to have the forcefield files correct so an MD simulation can be
generated with very little changes to the topol file after running
pdb2gmx. I do have a couple of questions that I hope you could help me
with. 

1) For grompp to compile successfully, must every molecule/residue have a
complete set of dihedral angle force constants?
2) It appears as though the residue to residue bond connectivity in my
force field files aren’t correctly parameterised. In NGLY-MOD-CGLY, a
bonded term is added to the topol for the C in NGLY to the N in CGLY. To
avoid this must I add an entry such as -C N in the [ bonds ] section of
the [ MOD1 ] residue definition in aminoacids.rtp? I’m unsure what “-C”
means.

Thanks for all your help.

Anthony


Dr Anthony Nash
Department of Chemistry
University College London





On 17/02/2016 21:46, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Mark Abraham" <gromacs.org_gmx-users-boun...@maillist.sys.kth.se
on behalf of mark.j.abra...@gmail.com> wrote:

>Hi,
>
>Your MOD residues are LYS connected by a ring, so the atom types and bonds
>sections should be pretty much as for LYS, with the ring parts perhaps
>along the lines of TYR. If you have done a full parameterization somehow
>and have specialized atom types, then yes, those bonded and non-bonded
>parameters need to go into the databases for grompp to look up (or then
>can
>also go in the [bonds] section of the .rtp, I think).
>
>Mark
>
>On Wed, Feb 17, 2016 at 10:39 PM Nash, Anthony <a.n...@ucl.ac.uk> wrote:
>
>>
>> Dear Mark,
>>
>>
>> I didn’t expect the problem was in ffnonbonded.itp. Problem solved.
>>Thanks
>> for the earlier hint.
>>
>> Anthony
>>
>> On 17/02/2016 18:01, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se
>>on
>> behalf of Nash, Anthony"
>> <gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
>> a.n...@ucl.ac.uk> wrote:
>>
>> >Hi Mark,
>> >
>> >Further to my earlier email. I’ve answered the query as to whether atom
>> >types and names can be the same. This is fine, judging by The N N and
>>O O
>> >in all of the amino acid types.
>> >
>> >I suspect I either have an additional problem or I have identified the
>> >root cause of my original problem. Every non-terminal amino acid has a
>>-C
>> >N line in [ bonds ]. My custom residue does not. I assume this is to do
>> >with amino acid connectivity. How do I implement this single bond if I
>> >have the connectivity NGLY-MOD-CGLY, where EVERY atom in MOD is
>>completely
>> >new? I assume this is where you meant use existing types. So in theory
>>I
>> >could just change the very first atom of MOD, which would be MN1 to N,
>>and
>> >then go into ffbonded and just change MN1-MC1 to N-MC1 for the first
>>bond
>> >parameter the MOD residue.
>> >
>> >Apologies for rambling on, I *think* I know what I am doing.
>> >
>> >Thanks
>> >Anthony
>> >
>> >Dr Anthony Nash
>> >Department of Chemistry
>> >University College London
>> >
>> >
>> >
>> >
>> >
>> >On 17/02/2016 17:18, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se
>> on
>> >behalf of Nash, Anthony"
>> ><gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
>> >a.n...@ucl.ac.uk> wrote:
>> >
>> >>Hi Mark,
>> >>
>> >>Thanks for the reply. I’m a little confused when you say “Choose
>>existing
>> >>types”. Are you saying that I am confined to only those atom types
>>that
>> >>come along with the forcefield upon installation, or that I can add
>>new
>> >>types but I’ve made a slight mishap between when specifying them
>>within
>> >>the [atom] section of my new residue?
>> >>
>> >>Also, I am guessing there is no issue with atom NAME and TYPE being
>> >>identical?
>> >>
>> >>Many thanks
>> >>Anthony
>> >>
>> >>Dr Anthony Nash
>> >>Department of Chemistry
>> >>University College London
>> >>
>> >>
>> >>
>> >>
>> >>
>> >>On 17/02/2016 16:43,
>>"gromacs.org_gmx-users-boun...@maillist.sys.kth.se
>> >>on
>> >>behalf of Mar

Re: [gmx-users] Unknown bond_atomtype

[Nash, Anthony]

Dear Mark,


I didn’t expect the problem was in ffnonbonded.itp. Problem solved. Thanks
for the earlier hint.

Anthony

On 17/02/2016 18:01, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Nash, Anthony"
<gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
a.n...@ucl.ac.uk> wrote:

>Hi Mark,
>
>Further to my earlier email. I’ve answered the query as to whether atom
>types and names can be the same. This is fine, judging by The N N and O O
>in all of the amino acid types.
>
>I suspect I either have an additional problem or I have identified the
>root cause of my original problem. Every non-terminal amino acid has a -C
>N line in [ bonds ]. My custom residue does not. I assume this is to do
>with amino acid connectivity. How do I implement this single bond if I
>have the connectivity NGLY-MOD-CGLY, where EVERY atom in MOD is completely
>new? I assume this is where you meant use existing types. So in theory I
>could just change the very first atom of MOD, which would be MN1 to N, and
>then go into ffbonded and just change MN1-MC1 to N-MC1 for the first bond
>parameter the MOD residue.
>
>Apologies for rambling on, I *think* I know what I am doing.
>
>Thanks
>Anthony
>
>Dr Anthony Nash
>Department of Chemistry
>University College London
>
>
>
>
>
>On 17/02/2016 17:18, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
>behalf of Nash, Anthony"
><gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
>a.n...@ucl.ac.uk> wrote:
>
>>Hi Mark,
>>
>>Thanks for the reply. I’m a little confused when you say “Choose existing
>>types”. Are you saying that I am confined to only those atom types that
>>come along with the forcefield upon installation, or that I can add new
>>types but I’ve made a slight mishap between when specifying them within
>>the [atom] section of my new residue?
>>
>>Also, I am guessing there is no issue with atom NAME and TYPE being
>>identical?
>>
>>Many thanks
>>Anthony
>>
>>Dr Anthony Nash
>>Department of Chemistry
>>University College London
>>
>>
>>
>>
>>
>>On 17/02/2016 16:43, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se
>>on
>>behalf of Mark Abraham"
>><gromacs.org_gmx-users-boun...@maillist.sys.kth.se
>>on behalf of mark.j.abra...@gmail.com> wrote:
>>
>>>Hi,
>>>
>>>You've specified a type for your atom in [atoms] and elsewhere a bond
>>>that
>>>uses it. Grompp has to find parameters for a bond between those two
>>>types,
>>>etc. Choose existing types ;-)
>>>
>>>Mark
>>>
>>>On Wed, 17 Feb 2016 17:27 Nash, Anthony <a.n...@ucl.ac.uk> wrote:
>>>
>>>> Hi all,
>>>>
>>>> As per a previous email (cross linking two peptide chains), I¹ve
>>>>created a
>>>> brand new crosslink (think disulphide bond) residue from scratch. I
>>>>have
>>>> defined it in all the files necessary (.rtp, residuetypes, specbond,
>>>> atomtypes, ffbonded, ffnonbonded) and it has got past pdb2gmx with no
>>>> problems.
>>>>
>>>> Unfortunately when I run grompp (5.0.4) I get the following error:
>>>>
>>>> Fatal error:
>>>> Unknown bond_atomtype MN1 <‹‹ first atom in my new residue
>>>>
>>>>
>>>> At first I thought my topology might be pointing to the wrong
>>>>forcefield,
>>>> but I¹ve checked and double checked:
>>>>
>>>> #include 
>>>>"/usr/local/gromacs/share/gromacs/top/amber99sb.ff/forcefield.itp"
>>>> #include "modic.itp"
>>>> ;#ifdef POSRES
>>>> ;#include "posre.itp"
>>>> ;#endif
>>>>
>>>> [ system ]
>>>> ; Name
>>>> MODIC with glycine terminal ends
>>>>
>>>> [ molecules ]
>>>> ; Compound#mols
>>>> MODIC 1
>>>>
>>>>
>>>> I¹m going to take a guess at what the problem is. Each of the atoms in
>>>>the
>>>> residue was derived from scratch.Therefore, they have a completely new
>>>> type. I wasn¹t feeling very creative with my naming convention so the
>>>>atom
>>>> name and the atom type are identical e.g.,
>>>>
>>>> ;MODIC crosslink
>>>> [ MOD1 ]
>>>>  [ atoms ]
>>>> ; NAMETYPECHARGE  NUMBER
>>>> MN1 MN1 -0.3640 

Re: [gmx-users] Unknown bond_atomtype

[Nash, Anthony]

Hi Mark,

Further to my earlier email. I’ve answered the query as to whether atom
types and names can be the same. This is fine, judging by The N N and O O
in all of the amino acid types.

I suspect I either have an additional problem or I have identified the
root cause of my original problem. Every non-terminal amino acid has a -C
N line in [ bonds ]. My custom residue does not. I assume this is to do
with amino acid connectivity. How do I implement this single bond if I
have the connectivity NGLY-MOD-CGLY, where EVERY atom in MOD is completely
new? I assume this is where you meant use existing types. So in theory I
could just change the very first atom of MOD, which would be MN1 to N, and
then go into ffbonded and just change MN1-MC1 to N-MC1 for the first bond
parameter the MOD residue.

Apologies for rambling on, I *think* I know what I am doing.

Thanks
Anthony

Dr Anthony Nash
Department of Chemistry
University College London





On 17/02/2016 17:18, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Nash, Anthony"
<gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
a.n...@ucl.ac.uk> wrote:

>Hi Mark,
>
>Thanks for the reply. I’m a little confused when you say “Choose existing
>types”. Are you saying that I am confined to only those atom types that
>come along with the forcefield upon installation, or that I can add new
>types but I’ve made a slight mishap between when specifying them within
>the [atom] section of my new residue?
>
>Also, I am guessing there is no issue with atom NAME and TYPE being
>identical?
>
>Many thanks
>Anthony
>
>Dr Anthony Nash
>Department of Chemistry
>University College London
>
>
>
>
>
>On 17/02/2016 16:43, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
>behalf of Mark Abraham" <gromacs.org_gmx-users-boun...@maillist.sys.kth.se
>on behalf of mark.j.abra...@gmail.com> wrote:
>
>>Hi,
>>
>>You've specified a type for your atom in [atoms] and elsewhere a bond
>>that
>>uses it. Grompp has to find parameters for a bond between those two
>>types,
>>etc. Choose existing types ;-)
>>
>>Mark
>>
>>On Wed, 17 Feb 2016 17:27 Nash, Anthony <a.n...@ucl.ac.uk> wrote:
>>
>>> Hi all,
>>>
>>> As per a previous email (cross linking two peptide chains), I¹ve
>>>created a
>>> brand new crosslink (think disulphide bond) residue from scratch. I
>>>have
>>> defined it in all the files necessary (.rtp, residuetypes, specbond,
>>> atomtypes, ffbonded, ffnonbonded) and it has got past pdb2gmx with no
>>> problems.
>>>
>>> Unfortunately when I run grompp (5.0.4) I get the following error:
>>>
>>> Fatal error:
>>> Unknown bond_atomtype MN1 <‹‹ first atom in my new residue
>>>
>>>
>>> At first I thought my topology might be pointing to the wrong
>>>forcefield,
>>> but I¹ve checked and double checked:
>>>
>>> #include 
>>>"/usr/local/gromacs/share/gromacs/top/amber99sb.ff/forcefield.itp"
>>> #include "modic.itp"
>>> ;#ifdef POSRES
>>> ;#include "posre.itp"
>>> ;#endif
>>>
>>> [ system ]
>>> ; Name
>>> MODIC with glycine terminal ends
>>>
>>> [ molecules ]
>>> ; Compound#mols
>>> MODIC 1
>>>
>>>
>>> I¹m going to take a guess at what the problem is. Each of the atoms in
>>>the
>>> residue was derived from scratch.Therefore, they have a completely new
>>> type. I wasn¹t feeling very creative with my naming convention so the
>>>atom
>>> name and the atom type are identical e.g.,
>>>
>>> ;MODIC crosslink
>>> [ MOD1 ]
>>>  [ atoms ]
>>> ; NAMETYPECHARGE  NUMBER
>>> MN1 MN1 -0.3640 1
>>> MC1 MC1 0.0358  2
>>> MC15MC150.4871  3
>>> Š
>>>
>>> Š
>>>
>>> Will this have confused grompp?
>>>
>>> Many thanks
>>> Anthony
>>>
>>> Dr Anthony Nash
>>> Department of Chemistry
>>> University College London
>>>
>>>
>>>
>>> --
>>> Gromacs Users mailing list
>>>
>>> * Please search the archive at
>>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>>> posting!
>>>
>>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>>
>>> * For (un)subscribe requests visit
>>> https://maillist.sys.kth.se/m

Re: [gmx-users] Unknown bond_atomtype

[Nash, Anthony]

Hi Mark,

Thanks for the reply. I’m a little confused when you say “Choose existing
types”. Are you saying that I am confined to only those atom types that
come along with the forcefield upon installation, or that I can add new
types but I’ve made a slight mishap between when specifying them within
the [atom] section of my new residue?

Also, I am guessing there is no issue with atom NAME and TYPE being
identical?

Many thanks
Anthony

Dr Anthony Nash
Department of Chemistry
University College London





On 17/02/2016 16:43, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Mark Abraham" <gromacs.org_gmx-users-boun...@maillist.sys.kth.se
on behalf of mark.j.abra...@gmail.com> wrote:

>Hi,
>
>You've specified a type for your atom in [atoms] and elsewhere a bond that
>uses it. Grompp has to find parameters for a bond between those two types,
>etc. Choose existing types ;-)
>
>Mark
>
>On Wed, 17 Feb 2016 17:27 Nash, Anthony <a.n...@ucl.ac.uk> wrote:
>
>> Hi all,
>>
>> As per a previous email (cross linking two peptide chains), I¹ve
>>created a
>> brand new crosslink (think disulphide bond) residue from scratch. I have
>> defined it in all the files necessary (.rtp, residuetypes, specbond,
>> atomtypes, ffbonded, ffnonbonded) and it has got past pdb2gmx with no
>> problems.
>>
>> Unfortunately when I run grompp (5.0.4) I get the following error:
>>
>> Fatal error:
>> Unknown bond_atomtype MN1 <‹‹ first atom in my new residue
>>
>>
>> At first I thought my topology might be pointing to the wrong
>>forcefield,
>> but I¹ve checked and double checked:
>>
>> #include 
>>"/usr/local/gromacs/share/gromacs/top/amber99sb.ff/forcefield.itp"
>> #include "modic.itp"
>> ;#ifdef POSRES
>> ;#include "posre.itp"
>> ;#endif
>>
>> [ system ]
>> ; Name
>> MODIC with glycine terminal ends
>>
>> [ molecules ]
>> ; Compound#mols
>> MODIC 1
>>
>>
>> I¹m going to take a guess at what the problem is. Each of the atoms in
>>the
>> residue was derived from scratch.Therefore, they have a completely new
>> type. I wasn¹t feeling very creative with my naming convention so the
>>atom
>> name and the atom type are identical e.g.,
>>
>> ;MODIC crosslink
>> [ MOD1 ]
>>  [ atoms ]
>> ; NAMETYPECHARGE  NUMBER
>> MN1 MN1 -0.3640 1
>> MC1 MC1 0.0358  2
>> MC15MC150.4871  3
>> Š
>>
>> Š
>>
>> Will this have confused grompp?
>>
>> Many thanks
>> Anthony
>>
>> Dr Anthony Nash
>> Department of Chemistry
>> University College London
>>
>>
>>
>> --
>> Gromacs Users mailing list
>>
>> * Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>> posting!
>>
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>> * For (un)subscribe requests visit
>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>> send a mail to gmx-users-requ...@gromacs.org.
>>
>-- 
>Gromacs Users mailing list
>
>* Please search the archive at
>http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>posting!
>
>* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
>* For (un)subscribe requests visit
>https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>send a mail to gmx-users-requ...@gromacs.org.

-- 
Gromacs Users mailing list

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

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

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

[gmx-users] Unknown bond_atomtype

[Nash, Anthony]

Hi all,

As per a previous email (cross linking two peptide chains), I¹ve created a
brand new crosslink (think disulphide bond) residue from scratch. I have
defined it in all the files necessary (.rtp, residuetypes, specbond,
atomtypes, ffbonded, ffnonbonded) and it has got past pdb2gmx with no
problems.  

Unfortunately when I run grompp (5.0.4) I get the following error:

Fatal error:
Unknown bond_atomtype MN1 <‹‹ first atom in my new residue


At first I thought my topology might be pointing to the wrong forcefield,
but I¹ve checked and double checked:

#include "/usr/local/gromacs/share/gromacs/top/amber99sb.ff/forcefield.itp"
#include "modic.itp"
;#ifdef POSRES
;#include "posre.itp"
;#endif

[ system ]
; Name
MODIC with glycine terminal ends

[ molecules ]
; Compound#mols
MODIC 1


I¹m going to take a guess at what the problem is. Each of the atoms in the
residue was derived from scratch.Therefore, they have a completely new
type. I wasn¹t feeling very creative with my naming convention so the atom
name and the atom type are identical e.g.,

;MODIC crosslink
[ MOD1 ]
 [ atoms ]
; NAMETYPECHARGE  NUMBER
MN1 MN1 -0.3640 1
MC1 MC1 0.0358  2
MC15MC150.4871  3
Š

Š

Will this have confused grompp?

Many thanks
Anthony

Dr Anthony Nash
Department of Chemistry
University College London



-- 
Gromacs Users mailing list

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

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

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

Re: [gmx-users] crosslinking polypeptide chains

[Nash, Anthony]

Thank Justin and Mark,

Apologies for not stripping out earlier content from my lazy “Reply”
email. It was a slight of hand.

It had crossed my mind to simply make two separate residues as you both
suggested. Although I was trying to make most of this interchangeable with
the Amber suite (this is an amber forcefield, but unlike my colleagues, I
prefer working in Gromacs), in which XLeap can do this.
However, with respects to Gromacs, what would happen to the charges? I
would essentially have two charge groups if I split them, both not an
integer, but when summed as a pair they would be 0 charge. Is this
acceptable?  

Many thanks
Anthony 

Dr Anthony Nash
Department of Chemistry
University College London





On 16/02/2016 16:11, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Mark Abraham" <gromacs.org_gmx-users-boun...@maillist.sys.kth.se
on behalf of mark.j.abra...@gmail.com> wrote:

>Hi,
>
>Please start new topics in new emails, rather than confusing the web
>archives with replies to digests :-)
>
>On Tue, Feb 16, 2016 at 5:00 PM Nash, Anthony <a.n...@ucl.ac.uk> wrote:
>
>> Hi all,
>>
>> When executing pdb2gmx I am getting a fatal error due to dangling
>>bonds. I
>> know that it will be down to how I¹ve organised the .pdb file, I¹m just
>> lacking in the experience with TERs, -chainsep and -merge to solve
>>this. I
>> would appreciate hints/tips/outright-solutions.
>>
>> My protein is very simple. It is a proof of concept for forcefield
>> parameters I derived myself from QM data (in house software I've
>> developed). There are two chains, both are the triplet NGLY-MOD-CGLY.
>>The
>> MOD residue is actually a lysine covalently bound to a carbon ring which
>> is then covalently bound to the neighbouring chain¹s lysine. It is being
>> treaded as one complete residue, MOD. E.g.,
>>
>> NGLY-MOD-CGLY
>>
>>   |
>> NGLY-MOD-CGLY
>>
>>
>> Just to reiterate, there is only one instance of the residue MOD, the
>> instance has two defined backbones (I.e., backbone as part of each
>>chain).
>
>
>It seems like arranging for your .pdb file to have chain separators that
>make gmx pdb2gmx -chainsep whatever -merge interactive would work
>smoothly,
>and then have specbond mechanism do the crosslink. What can't work is a
>single MOD residue, because the early stages of pdb2gmx only know how to
>recognize (and handle termini for) an unbranched chain. You might need a
>(complete) LYS+ring residue in one chain, and a normal LYS in the other.
>
>Mark
>-- 
>Gromacs Users mailing list
>
>* Please search the archive at
>http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>posting!
>
>* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
>* For (un)subscribe requests visit
>https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>send a mail to gmx-users-requ...@gromacs.org.

-- 
Gromacs Users mailing list

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

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

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

[gmx-users] crosslinking polypeptide chains

[Nash, Anthony]

Hi all,

When executing pdb2gmx I am getting a fatal error due to dangling bonds. I
know that it will be down to how I¹ve organised the .pdb file, I¹m just
lacking in the experience with TERs, -chainsep and -merge to solve this. I
would appreciate hints/tips/outright-solutions.

My protein is very simple. It is a proof of concept for forcefield
parameters I derived myself from QM data (in house software I've
developed). There are two chains, both are the triplet NGLY-MOD-CGLY. The
MOD residue is actually a lysine covalently bound to a carbon ring which
is then covalently bound to the neighbouring chain¹s lysine. It is being
treaded as one complete residue, MOD. E.g.,

NGLY-MOD-CGLY

  |
NGLY-MOD-CGLY


Just to reiterate, there is only one instance of the residue MOD, the
instance has two defined backbones (I.e., backbone as part of each chain).

I would appreciate any assistance. Many thanks
Anthony



Dr Anthony Nash
Department of Chemistry
University College London





On 16/02/2016 15:46, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Mark Abraham"  wrote:

>Hi,
>
>On Tue, Feb 16, 2016 at 4:07 PM Nicolas Cheron <
>nicolas.cheron.bou...@gmail.com> wrote:
>
>> Dear all,
>>
>> Is there a place where the planned/hoped new features of Gromacs are
>> listed?
>
>
>Kind of. Many times people do some preliminary discussion on Redmine, and
>there's a large amount of stuff people have wished for e.g. at
>http://redmine.gromacs.org/projects/gromacs/roadmap, but also quite a few
>projects around that don't currently have a place there. We do releases
>because it's that time of year, rather than because a feature is thought
>to
>be ready, so there's not a natural list of things people are working on
>for
>that release.
>
>
>> I am wondering if there is a chance that EVB, HREX (without the
>> need to use plumed) or new QM/MM interface will be available at some
>>point.
>>
>
>I don't know of anybody with short-term plans for these - getting funding
>for software development, and people suited to do it, is always
>challenging, and so a lot of what is available gets channeled into core
>infrastructure, because that pays off for the whole user community. That
>said, there are some forms of HREX you can already do with standard
>GROMACS, and there is infrastructure work going on that would eventually
>permit tools like PLUMED to be implemented in more efficient ways. I know
>that the CPMD team is working on replacing their interface to GROMOS with
>one that will use GROMACS; for that, CPMD will do the driving. The QM-MM
>interface in GROMACS seems not to be interesting for anybody with time to
>work on it :-(.
>
>Mark
>
>I know that the developers are busy, but I am wondering what is in their
>> mind.
>>
>> Thank you.
>>
>> Nicolas
>> --
>> Gromacs Users mailing list
>>
>> * Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>> posting!
>>
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>> * For (un)subscribe requests visit
>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>> send a mail to gmx-users-requ...@gromacs.org.
>>
>-- 
>Gromacs Users mailing list
>
>* Please search the archive at
>http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>posting!
>
>* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
>* For (un)subscribe requests visit
>https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>send a mail to gmx-users-requ...@gromacs.org.

-- 
Gromacs Users mailing list

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

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

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

Re: [gmx-users] membed in mdrun VERSION 5.0.4

[Nash, Anthony]

My apologies that this hasn’t made it into one large message. I would like
to add two further points. It appears as though mdrun (serial) does not
output a .cpt file when -membed is specified. I just ran out of clock
time, and the check point file is just not there.


Also, in response to your earlier comment about using gmx convert-tpr,
-membed appears to remove some lipids and solvent molecules but it won’t
tell you which ones until you have the final .gro file. Therefore, I don’t
believe I can make a .ndx to use with gmx-convert-tpr.

Thanks
Anthony  

On 09/01/2016 09:28, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Nash, Anthony"
<gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
a.n...@ucl.ac.uk> wrote:

>Hi,
>Just to add to my earlier message, I went through all release notes after
>5.0.4, and besides a change in membed documentation, I can’t see a
>reference to a parallelized version of mdrun -membed.
>
>
>Thanks
>Anthony
>
>On 08/01/2016 23:36, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
>behalf of Nash, Anthony"
><gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
>a.n...@ucl.ac.uk> wrote:
>
>>Justin and Mark, many thanks for your help.
>>
>>With regards to the parallelization, when did parallel membed become
>>supported? I¹ve just tried on 5.0.4 and I¹m getting the response:
>>
>>Reading file membed_NPT_B.tpr, VERSION 5.0.4 (double precision)
>>
>>---
>>Program mdrun_mpi_d, VERSION 5.0.4
>>Source code file:
>>/dev/shm/tmp.3LEFQsSzrx/gromacs-5.0.4/src/programs/mdrun/membed.c, line:
>>1056
>>
>>
>>Input error or input inconsistency:
>>Sorry, parallel g_membed is not yet fully functional.
>>For more information and tips for troubleshooting, please check the
>>GROMACS
>>website at http://www.gromacs.org/Documentation/Errors
>>---
>>
>>Error on rank 0, will try to stop all ranks
>>Halting parallel program mdrun_mpi_d on CPU 0 out of 2
>>
>>
>>
>>Oh course, I understand this could be answered by installing the latest
>>version, but I¹m ruling out anything I¹ve done wrong before asking
>>central
>>IT to put up a test module.
>>
>>Thanks for the help
>>Anthony 
>>
>>
>>
>>
>>
>>On 08/01/2016 13:24, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se
>>on
>>behalf of Mark Abraham"
>><gromacs.org_gmx-users-boun...@maillist.sys.kth.se
>>on behalf of mark.j.abra...@gmail.com> wrote:
>>
>>>Hi,
>>>
>>>The current implementation of membed is quite sane and should support
>>>all
>>>kinds of parallelization. (It is difficult to say much that is nice
>>>about
>>>the first implementation, though!)
>>>
>>>gmx trjconv and convert-tpr are good for making subsets from well chosen
>>>index groups.
>>>
>>>Mark
>>>
>>>On Fri, 8 Jan 2016 13:35 Justin Lemkul <jalem...@vt.edu> wrote:
>>>
>>>>
>>>>
>>>> On 1/8/16 3:38 AM, Nash, Anthony wrote:
>>>> > Many thanks Justin, that¹s solved it.
>>>> >
>>>> > The simulation is now running, reporting that 8 POPC molecules and
>>>>12
>>>>SOL
>>>> > molecules have been removed. However, I have a bit of a problem. I
>>>>don¹t
>>>> > seem to be able to get the -membed option working on a parallel
>>>> > installation of gromacs, only serial. I assume parallel isn¹t
>>>>supported?
>>>> >
>>>>
>>>> No idea, but the -membed incorporation into mdrun is currently a very
>>>>big
>>>> hack
>>>> to make it work, so you should expect limited support for features.
>>>>
>>>> > Secondly, with my .tpr file very different from my running .trr (due
>>>>to
>>>> > removed molecules), I don¹t know of a way of pulling out a single
>>>>frame
>>>> > from the trajectory to see how the simulation is progressing. (I
>>>>would
>>>> > have just taken the first frame as a .gro, downloaded the latest
>>>>.trr,
>>>> and
>>>> > thrown them both into VMD). Any thoughts?
>>>> >
>>>>
>>>> Does the program report which residues were removed?  If so, you can
>>>>just
>>>> remove
>>>> those from the coord

Re: [gmx-users] membed in mdrun VERSION 5.0.4

[Nash, Anthony]

Hi,
Just to add to my earlier message, I went through all release notes after
5.0.4, and besides a change in membed documentation, I can’t see a
reference to a parallelized version of mdrun -membed.


Thanks
Anthony

On 08/01/2016 23:36, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Nash, Anthony"
<gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
a.n...@ucl.ac.uk> wrote:

>Justin and Mark, many thanks for your help.
>
>With regards to the parallelization, when did parallel membed become
>supported? I¹ve just tried on 5.0.4 and I¹m getting the response:
>
>Reading file membed_NPT_B.tpr, VERSION 5.0.4 (double precision)
>
>---
>Program mdrun_mpi_d, VERSION 5.0.4
>Source code file: 
>/dev/shm/tmp.3LEFQsSzrx/gromacs-5.0.4/src/programs/mdrun/membed.c, line:
>1056
>
>
>Input error or input inconsistency:
>Sorry, parallel g_membed is not yet fully functional.
>For more information and tips for troubleshooting, please check the
>GROMACS
>website at http://www.gromacs.org/Documentation/Errors
>---
>
>Error on rank 0, will try to stop all ranks
>Halting parallel program mdrun_mpi_d on CPU 0 out of 2
>
>
>
>Oh course, I understand this could be answered by installing the latest
>version, but I¹m ruling out anything I¹ve done wrong before asking central
>IT to put up a test module.
>
>Thanks for the help
>Anthony 
>
>
>
>
>
>On 08/01/2016 13:24, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
>behalf of Mark Abraham" <gromacs.org_gmx-users-boun...@maillist.sys.kth.se
>on behalf of mark.j.abra...@gmail.com> wrote:
>
>>Hi,
>>
>>The current implementation of membed is quite sane and should support all
>>kinds of parallelization. (It is difficult to say much that is nice about
>>the first implementation, though!)
>>
>>gmx trjconv and convert-tpr are good for making subsets from well chosen
>>index groups.
>>
>>Mark
>>
>>On Fri, 8 Jan 2016 13:35 Justin Lemkul <jalem...@vt.edu> wrote:
>>
>>>
>>>
>>> On 1/8/16 3:38 AM, Nash, Anthony wrote:
>>> > Many thanks Justin, that¹s solved it.
>>> >
>>> > The simulation is now running, reporting that 8 POPC molecules and 12
>>>SOL
>>> > molecules have been removed. However, I have a bit of a problem. I
>>>don¹t
>>> > seem to be able to get the -membed option working on a parallel
>>> > installation of gromacs, only serial. I assume parallel isn¹t
>>>supported?
>>> >
>>>
>>> No idea, but the -membed incorporation into mdrun is currently a very
>>>big
>>> hack
>>> to make it work, so you should expect limited support for features.
>>>
>>> > Secondly, with my .tpr file very different from my running .trr (due
>>>to
>>> > removed molecules), I don¹t know of a way of pulling out a single
>>>frame
>>> > from the trajectory to see how the simulation is progressing. (I
>>>would
>>> > have just taken the first frame as a .gro, downloaded the latest
>>>.trr,
>>> and
>>> > thrown them both into VMD). Any thoughts?
>>> >
>>>
>>> Does the program report which residues were removed?  If so, you can
>>>just
>>> remove
>>> those from the coordinate file you had before.  Otherwise, no, probably
>>> not.
>>> Like I said, big hack to get it working in the present version.
>>>
>>> -Justin
>>>
>>> --
>>> ==
>>>
>>> Justin A. Lemkul, Ph.D.
>>> Ruth L. Kirschstein NRSA Postdoctoral Fellow
>>>
>>> Department of Pharmaceutical Sciences
>>> School of Pharmacy
>>> Health Sciences Facility II, Room 629
>>> University of Maryland, Baltimore
>>> 20 Penn St.
>>> Baltimore, MD 21201
>>>
>>> jalem...@outerbanks.umaryland.edu | (410) 706-7441
>>> http://mackerell.umaryland.edu/~jalemkul
>>>
>>> ==
>>> --
>>> Gromacs Users mailing list
>>>
>>> * Please search the archive at
>>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>>> posting!
>>>
>>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>>
>>> * For (un)subscribe requests visit
>>> https://maillist.sys.kth.se/mailman/list

Re: [gmx-users] membed in mdrun VERSION 5.0.4

[Nash, Anthony]

Many thanks Justin, that’s solved it.

The simulation is now running, reporting that 8 POPC molecules and 12 SOL
molecules have been removed. However, I have a bit of a problem. I don’t
seem to be able to get the -membed option working on a parallel
installation of gromacs, only serial. I assume parallel isn’t supported?

Secondly, with my .tpr file very different from my running .trr (due to
removed molecules), I don’t know of a way of pulling out a single frame
from the trajectory to see how the simulation is progressing. (I would
have just taken the first frame as a .gro, downloaded the latest .trr, and
thrown them both into VMD). Any thoughts?

Many thanks
Anthony 




On 06/01/2016 19:18, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Justin Lemkul"
<gromacs.org_gmx-users-boun...@maillist.sys.kth.se on behalf of
jalem...@vt.edu> wrote:

>
>
>On 1/6/16 2:00 PM, Nash, Anthony wrote:
>> Hi all,
>>
>> I thought I would try using the -membed option of mdrun to insert a
>> helical dimer into a lipid bilayer. I¹ve followed the .mdp instructions
>>on
>> g_membed to generate my required .tpr file. Upon calling grommp I get:
>>
>> ERROR 1 [file membed_NPT.mdp]:
>>Energy group exclusions are not (yet) implemented for the Verlet
>>scheme
>>
>>
>> My input file can be seen below. If seems as though the Integrator is
>>the
>> issue, even though I¹ve picked the correct one. Is a place holder
>>feature?
>>
>
>It's not the integrator that's the problem.  It is, as the error states,
>the use 
>of the Verlet cutoff scheme.  Presumably "cutoff-scheme = group" would
>circumvent this issue.
>
>-Justin
>
>> Many thanks
>> Anthony
>>
>> integrator = md
>> energygrps = Protein
>> freezegrps = Protein
>> freezedim = Y Y Y
>> energygrp_excl = Protein Protein
>>
>>
>> nsteps  = 500   ; 2 * 50 = 1000 ps (2 ns)
>> dt  = 0.002 ; 2 fs
>> nstxout = 1 ; save coordinates every 0.2 ps
>> nstvout = 1 ; save velocities every 0.2 ps
>> nstenergy   = 1 ; save energies every 0.2 ps
>> nstlog  = 1 ; update log file every 0.2 ps
>>
>> continuation= yes   ; Restarting after NVT
>> constraint_algorithm = lincs; holonomic constraints
>> constraints = all-bonds ; all bonds (even heavy atom-H
>> bonds) constrained
>> lincs_iter  = 1 ; accuracy of LINCS
>> lincs_order = 4 ; also related to accuracy
>>
>> ns_type = grid  ; search neighboring grid cels
>> nstlist = 5 ; 10 fs
>> rlist   = 1.0   ; short-range neighborlist cutoff (in
>>nm)
>> rcoulomb= 1.0   ; short-range electrostatic cutoff (in
>>nm)
>> rvdw= 1.0   ; short-range van der Waals cutoff (in
>>nm)
>>
>> coulombtype = PME   ; Particle Mesh Ewald for long-range
>> electrostatics
>> pme_order   = 4 ; cubic interpolation
>> fourierspacing  = 0.16  ; grid spacing for FFT
>>
>> tcoupl  = Berendsen   ;Nose-Hoover  ; More
>> accurate thermostat
>> tc-grps = Protein  POPC SOL ; three coupling groups - more
>> accurate
>> tau_t   = 0.5   0.5 0.5 ; time constant, in ps
>> ref_t   = 310   310 310 ; reference temperature,
>> one for each group, in K
>>
>> pcoupl  = Parrinello-Rahman ; Pressure coupling on in
>>NPT
>> pcoupltype  = semiisotropic ; uniform scaling of x-y box
>> vectors, independent z
>> tau_p   = 5.0   ; time constant, in ps
>> ref_p   = 1.0   1.0 ; reference pressure,
>>x-y,
>> z (in bar)
>> compressibility = 4.5e-54.5e-5  ; isothermal compressibility,
>> bar^-1
>>
>> pbc = xyz   ; 3-D PBC
>> DispCorr= EnerPres  ; account for cut-off vdW scheme
>> gen_vel = no; Velocity generation is off
>> nstcomm = 1
>> comm-mode   = Linear
>> comm-grps   = POPC_Protein SOL
>> refcoord_scaling = com
>>
>>
>
>-- 
>==
>
>Justin A. Lemkul, Ph.D.
>Ruth L. Kirschstein NRSA Postdoctoral Fellow
>
>Department of Pharmaceutical Sciences
>School of Pharmacy
>Health Sci

Re: [gmx-users] membed in mdrun VERSION 5.0.4

[Nash, Anthony]

Justin and Mark, many thanks for your help.

With regards to the parallelization, when did parallel membed become
supported? I¹ve just tried on 5.0.4 and I¹m getting the response:

Reading file membed_NPT_B.tpr, VERSION 5.0.4 (double precision)

---
Program mdrun_mpi_d, VERSION 5.0.4
Source code file: 
/dev/shm/tmp.3LEFQsSzrx/gromacs-5.0.4/src/programs/mdrun/membed.c, line:
1056


Input error or input inconsistency:
Sorry, parallel g_membed is not yet fully functional.
For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors
---

Error on rank 0, will try to stop all ranks
Halting parallel program mdrun_mpi_d on CPU 0 out of 2



Oh course, I understand this could be answered by installing the latest
version, but I¹m ruling out anything I¹ve done wrong before asking central
IT to put up a test module.

Thanks for the help
Anthony 





On 08/01/2016 13:24, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Mark Abraham" <gromacs.org_gmx-users-boun...@maillist.sys.kth.se
on behalf of mark.j.abra...@gmail.com> wrote:

>Hi,
>
>The current implementation of membed is quite sane and should support all
>kinds of parallelization. (It is difficult to say much that is nice about
>the first implementation, though!)
>
>gmx trjconv and convert-tpr are good for making subsets from well chosen
>index groups.
>
>Mark
>
>On Fri, 8 Jan 2016 13:35 Justin Lemkul <jalem...@vt.edu> wrote:
>
>>
>>
>> On 1/8/16 3:38 AM, Nash, Anthony wrote:
>> > Many thanks Justin, that¹s solved it.
>> >
>> > The simulation is now running, reporting that 8 POPC molecules and 12
>>SOL
>> > molecules have been removed. However, I have a bit of a problem. I
>>don¹t
>> > seem to be able to get the -membed option working on a parallel
>> > installation of gromacs, only serial. I assume parallel isn¹t
>>supported?
>> >
>>
>> No idea, but the -membed incorporation into mdrun is currently a very
>>big
>> hack
>> to make it work, so you should expect limited support for features.
>>
>> > Secondly, with my .tpr file very different from my running .trr (due
>>to
>> > removed molecules), I don¹t know of a way of pulling out a single
>>frame
>> > from the trajectory to see how the simulation is progressing. (I would
>> > have just taken the first frame as a .gro, downloaded the latest .trr,
>> and
>> > thrown them both into VMD). Any thoughts?
>> >
>>
>> Does the program report which residues were removed?  If so, you can
>>just
>> remove
>> those from the coordinate file you had before.  Otherwise, no, probably
>> not.
>> Like I said, big hack to get it working in the present version.
>>
>> -Justin
>>
>> --
>> ==
>>
>> Justin A. Lemkul, Ph.D.
>> Ruth L. Kirschstein NRSA Postdoctoral Fellow
>>
>> Department of Pharmaceutical Sciences
>> School of Pharmacy
>> Health Sciences Facility II, Room 629
>> University of Maryland, Baltimore
>> 20 Penn St.
>> Baltimore, MD 21201
>>
>> jalem...@outerbanks.umaryland.edu | (410) 706-7441
>> http://mackerell.umaryland.edu/~jalemkul
>>
>> ==
>> --
>> Gromacs Users mailing list
>>
>> * Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>> posting!
>>
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>> * For (un)subscribe requests visit
>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>> send a mail to gmx-users-requ...@gromacs.org.
>-- 
>Gromacs Users mailing list
>
>* Please search the archive at
>http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>posting!
>
>* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
>* For (un)subscribe requests visit
>https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>send a mail to gmx-users-requ...@gromacs.org.

-- 
Gromacs Users mailing list

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

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

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

[gmx-users] membed in mdrun VERSION 5.0.4

[Nash, Anthony]

Hi all,

I thought I would try using the -membed option of mdrun to insert a
helical dimer into a lipid bilayer. I¹ve followed the .mdp instructions on
g_membed to generate my required .tpr file. Upon calling grommp I get:

ERROR 1 [file membed_NPT.mdp]:
  Energy group exclusions are not (yet) implemented for the Verlet scheme


My input file can be seen below. If seems as though the Integrator is the
issue, even though I¹ve picked the correct one. Is a place holder feature?

Many thanks
Anthony

integrator = md
energygrps = Protein
freezegrps = Protein
freezedim = Y Y Y
energygrp_excl = Protein Protein


nsteps  = 500   ; 2 * 50 = 1000 ps (2 ns)
dt  = 0.002 ; 2 fs
nstxout = 1 ; save coordinates every 0.2 ps
nstvout = 1 ; save velocities every 0.2 ps
nstenergy   = 1 ; save energies every 0.2 ps
nstlog  = 1 ; update log file every 0.2 ps

continuation= yes   ; Restarting after NVT
constraint_algorithm = lincs; holonomic constraints
constraints = all-bonds ; all bonds (even heavy atom-H
bonds) constrained
lincs_iter  = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy

ns_type = grid  ; search neighboring grid cels
nstlist = 5 ; 10 fs
rlist   = 1.0   ; short-range neighborlist cutoff (in nm)
rcoulomb= 1.0   ; short-range electrostatic cutoff (in nm)
rvdw= 1.0   ; short-range van der Waals cutoff (in nm)

coulombtype = PME   ; Particle Mesh Ewald for long-range
electrostatics
pme_order   = 4 ; cubic interpolation
fourierspacing  = 0.16  ; grid spacing for FFT

tcoupl  = Berendsen   ;Nose-Hoover  ; More
accurate thermostat
tc-grps = Protein  POPC SOL ; three coupling groups - more
accurate
tau_t   = 0.5   0.5 0.5 ; time constant, in ps
ref_t   = 310   310 310 ; reference temperature,
one for each group, in K

pcoupl  = Parrinello-Rahman ; Pressure coupling on in NPT
pcoupltype  = semiisotropic ; uniform scaling of x-y box
vectors, independent z
tau_p   = 5.0   ; time constant, in ps
ref_p   = 1.0   1.0 ; reference pressure, x-y,
z (in bar)
compressibility = 4.5e-54.5e-5  ; isothermal compressibility,
bar^-1

pbc = xyz   ; 3-D PBC
DispCorr= EnerPres  ; account for cut-off vdW scheme
gen_vel = no; Velocity generation is off
nstcomm = 1
comm-mode   = Linear
comm-grps   = POPC_Protein SOL
refcoord_scaling = com


-- 
Gromacs Users mailing list

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

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

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

[gmx-users] inserting TM protein dimer into lipid bilayer using Gromacs

[Nash, Anthony]

Dear all,

It¹s been almost two and a 1/2 years since I tried my hands at TM protein
modelling using Gromacs. What is the latest and most reliable means of
inserting a TM alpha helical dimer into a lipid bilayer using Grimaces (if
possible)?

Thanks
Anthony 

Dr Anthony Nash
Department of Chemistry
University College London


>

-- 
Gromacs Users mailing list

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

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

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

Re: [gmx-users] gmx mdrun std::bad_alloc whilst using PLUMED

[Nash, Anthony]

Thanks Mark,

I threw an email across to the plumed group this morning. I was surprised
to get a reply almost immediately. It *could* be the memory allocation
required to define the grid spacing in PLUMED.

Thanks
Anthony

Dr Anthony Nash
Department of Chemistry
University College London





On 17/11/2015 22:24, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Mark Abraham" <gromacs.org_gmx-users-boun...@maillist.sys.kth.se
on behalf of mark.j.abra...@gmail.com> wrote:

>Hi,
>
>GROMACS is apparently the first to notice that memory is a problem, but
>you
>should also be directing questions about memory use with different kinds
>of
>CVs to the PLUMED people. mdrun knows nothing at all about the PLUMED CVs.
>The most likely explanation is that they have some data structure that
>works OK on small scale problems, but which doesn't do well as the number
>of atoms, CVs, CV complexity, and/or ranks increases.
>
>Mark
>
>On Tue, Nov 17, 2015 at 11:05 PM Nash, Anthony <a.n...@ucl.ac.uk> wrote:
>
>> Hi all,
>>
>> I am using PLUMED 2.2 and gromacs 5.0.4. For a while I had been testing
>> the viability of three collective variables for plumed, two dihedral
>> angles and one centre of mass distance. After observing my dimer rotate
>> about each other I decided it needed an intrahelical distance between
>>two
>> of the dihedral atoms (A,B,C,D), per helix, to sample the CV space
>>whilst
>> maintaining the Œregular¹ alpha-helical structure (the dihedral sampling
>> was coming from the protein uncoiling rather than rotating). Note: it is
>> likely that I will change these distances to the built-in alpha helical
>> CV.
>>
>> The moment I increased the number of CVs from three to five, gromacs
>> throws out a memory error. When I remove the 5th CV it still crashes.
>>When
>> I remove the 4th it stops crashing.
>>
>> ‹‹‹
>> CLUSTER OUTPUT FILE
>> ‹‹‹
>>
>>
>> starting mdrun 'NEU_MUT in POPC in water'
>> 5000 steps, 10.0 ps.
>>
>> ---
>> Program: gmx mdrun, VERSION 5.0.4
>>
>> Memory allocation failed:
>> std::bad_alloc
>>
>> For more information and tips for troubleshooting, please check the
>>GROMACS
>> website at http://www.gromacs.org/Documentation/Errors
>> ---
>> Halting parallel program mdrun_mpi_d on CPU 0 out of 12
>>
>>
>>
>>
>> It halts all 12 processes and the job dies. I increased the memory so I
>>am
>> using 43.2 GB of RAM distributed over 12 processes. The job still fails
>> (but then, allocation of memory is very different to not having any
>>memory
>> at all).
>>
>> The contents of the gromacs.log file only reports the initialisation of
>> gromacs followed by the initialisation of plumed. After this I would
>>have
>> expected the regular MD stepping output. I¹ve included the plumed
>> initialisation below. I would appreciate any help. I suspect the problem
>> lies with the 4th and 5th CV although systematically removing them and
>> playing around with the parameters hasn¹t yielded anything yet. Please
>> ignore what parameter values I have set. Besides the atom number
>> everything else is a result of me trying to work out where certain
>>ranges
>> of values is causing PLUMED to exit and gromacs to crash. PLUMED input
>> file below:
>>
>>
>> ‹‹‹
>> PLUMED INPUTFILE
>> ‹‹‹
>>
>> phi: TORSION ATOMS=214,230,938,922
>> psi: TORSION ATOMS=785,801,367,351
>>
>> c1: COM ATOMS=1-571
>> c2: COM ATOMS=572-1142
>> COMdist: DISTANCE ATOMS=c1,c2
>>
>> d1: DISTANCE ATOMS=214,367
>> d2: DISTANCE ATOMS=938,785
>>
>> UPPER_WALLS ARG=COMdist AT=2.5 KAPPA=1000 EXP=2.0 EPS=1.0 OFFSET=0
>> LABEL=COMuwall
>> LOWER_WALLS ARG=COMdist AT=1.38 KAPPA=1000 EXP=2.0 EPS=1.0 OFFSET=0
>> LABEL=COMlwall
>>
>> UPPER_WALLS ARG=d1 AT=1.260 KAPPA=1000 EXP=2.0 EPS=1.0 OFFSET=0
>> LABEL=d1uwall
>> LOWER_WALLS ARG=d1 AT=1.228 KAPPA=1000 EXP=2.0 EPS=1.0 OFFSET=0
>> LABEL=d1lwall
>>
>> UPPER_WALLS ARG=d2 AT=1.228 KAPPA=1000 EXP=2.0 EPS=1.0 OFFSET=0
>> LABEL=d2uwall
>> LOWER_WALLS ARG=d2 AT=1.196 KAPPA=1000 EXP=2.0 EPS=1.0 OFFSET=0
>> LABEL=d2lwall
>>
>> METAD ...
>> LABEL=metad
>> ARG=phi,psi,COMdist,d1,d2
>> PACE=1
>> HEIGHT=0.2
>> SIGMA=0.06,0.06,0.06,0.06,0.06
>> FILE=HIL

Re: [gmx-users] restart issues with Gromacs

[Nash, Anthony]

Thanks Mark,

Ok, so it sounds very much like it¹s on the cluster side. I¹ll fire this
across to one of the sys admins and see if they can find out what the
problem is, although I have no idea what and if they will find anything.
>From the looks of it, the reading of the checkpoint file data was fine
(else I expect an MD5 hashkey error IF my memory serves me right).

Thanks again
Anthony

Dr Anthony Nash
Department of Chemistry
University College London





On 15/11/2015 22:12, "gromacs.org_gmx-users-boun...@maillist.sys.kth.se on
behalf of Mark Abraham" <gromacs.org_gmx-users-boun...@maillist.sys.kth.se
on behalf of mark.j.abra...@gmail.com> wrote:

>Hi,
>
>GROMACS assumes the file systems actually write to disk when you use the
>system call that means that, and works correctly if so. But if the file
>system or its configuration don't actually do that (for "performance" or
>erroneous reasons), then all bets are off. mdrun can't even know if it's
>being lied to, because, well, it's being lied to...
>
>Mark
>
>On Sun, 15 Nov 2015 22:30 Nash, Anthony <a.n...@ucl.ac.uk> wrote:
>
>> Hi all,
>>
>> Running Grimaces 5.0.4 with PLUMED 2.2 on a cluster, number of ranks
>>(MPI
>> processes) is 24. The simulation successfully ran for the maximum
>>cluster
>> wall time (48 hours).
>>
>> I attempt to restart the simulations using the following command (with a
>> sun microsystem grid engine submission script):
>>
>> gerun mdrun_mpi_d -deffnm neu_mut_meta_K -cpi neu_mut_meta_K.cpt
>>-noappend
>> -plumed
>>
>>
>> However, whilst the queue is telling me that the job is running, the
>> *.part0002.log file seems stuck at:
>>
>> ---
>> When dynamic load balancing gets turned on, these settings will change
>>to:
>> The maximum number of communication pulses is: X 1 Y 1
>> The minimum size for domain decomposition cells is 1.025 nm
>> The requested allowed shrink of DD cells (option -dds) is: 0.80
>> The allowed shrink of domain decomposition cells is: X 0.43 Y 0.56
>> The maximum allowed distance for charge groups involved in interactions
>>is:
>>  non-bonded interactions   1.025 nm
>> two-body bonded interactions  (-rdd)   1.025 nm
>>   multi-body bonded interactions  (-rdd)   1.025 nm
>>   atoms separated by up to 5 constraints  (-rcon)  1.025 nm
>> 
>>
>>
>> The cluster error file and output file (not gromacs file) contains no
>> warnings or errors. The gromacs log file contains no warnings or errors.
>>
>> I have seen this behaviour quite a number of times, going back to early
>> versions of gromacs 4 around late 2010 (I think). I got into the habit
>>of
>> a) copying backing up the .cpt files, and b) always using -noappend
>>option
>> to preserve the .trr file. Has there ever been an explanation as to why
>> this is happening?
>>
>> Many thanks
>> Anthony
>>
>> --
>> Gromacs Users mailing list
>>
>> * Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>> posting!
>>
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>> * For (un)subscribe requests visit
>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>> send a mail to gmx-users-requ...@gromacs.org.
>>
>-- 
>Gromacs Users mailing list
>
>* Please search the archive at
>http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>posting!
>
>* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
>* For (un)subscribe requests visit
>https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>send a mail to gmx-users-requ...@gromacs.org.

-- 
Gromacs Users mailing list

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

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

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

Re: [gmx-users] syntax for gmx distance

[Nash, Anthony]

Thanks Teemu,

That was spot on. I simply took the -select argument and dropped it into a
text file then used the -sf option instead. Works perfectly.

Thanks again
Anthony 




On 14/09/2015 12:34, "Nash, Anthony" <a.n...@ucl.ac.uk> wrote:

>Thanks Teemu,
>
>I’ll look into this further.
>
>Thanks
>Anthony 
>
>
>
>On 14/09/2015 12:32, "Teemu Murtola" <teemu.murt...@gmail.com> wrote:
>
>>The error message really looks like that your shell does not understand
>>the
>>quoting. You need to pass the whole selection as a single command-line
>>parameter (which your quotes in principle would do), but the error
>>message
>>indicates that each word is getting parsed as a separate selection, so
>>your
>>shell is passing them as separate arguments.
>>
>>Best regards,
>>Teemu
>>
>>On Mon, Sep 14, 2015, 11:28 Nash, Anthony <a.n...@ucl.ac.uk> wrote:
>>
>>Dear all,
>>
>>I understand that this is quite a basic question, but I think I need a
>>fresh set of eyes to figure out what¹s going on with gmx distance -select
>>syntax.
>>
>>In my index file, I have 23 groups, two of which are ³ZINC² and
>>³CARBONYL²
>>
>>Based on numerous mail archive suggestions and the gromacs website change
>>in tools page, I¹m using:
>>
>>$CHEM/rungromacs505 gmx distance -n index.ndx -f umb_1_all.trr -s
>>umb_1_umb.tpr -oh umb_1_distance_hist.xvg -oall umb_1_distance.xvg
>>-select
>>'com of group "ZINC" plus com of group "CARBONYL"'
>>
>>Which kicks up the error:
>>
>>Invalid command-line options
>>  In command-line option -select
>>In keyword 'com'
>>  required parameter 'of' not specified
>>invalid selection 'com'
>>  In command-line option -select
>>syntax error
>>invalid selection 'group'
>>  In command-line option -select
>>syntax error
>>invalid selection 'plus'
>>  In command-line option -select
>>In keyword 'com'
>>  required parameter 'of' not specified
>>invalid selection 'com'
>>  In command-line option -select
>>syntax error
>>invalid selection 'group'
>>
>>I¹d appreciate any suggestions. I¹ve also tried to do *something* via the
>>interactive command line (when -select is not supplied), but what it¹s
>>asking for is a bit cryptic.
>>
>>Thanks
>>Anthony
>>-- 
>>Gromacs Users mailing list
>>
>>* Please search the archive at
>>http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>>posting!
>>
>>* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>>* For (un)subscribe requests visit
>>https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>>send a mail to gmx-users-requ...@gromacs.org.
>
>-- 
>Gromacs Users mailing list
>
>* Please search the archive at
>http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>posting!
>
>* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
>* For (un)subscribe requests visit
>https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>send a mail to gmx-users-requ...@gromacs.org.

-- 
Gromacs Users mailing list

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

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

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

[gmx-users] syntax for gmx distance

[Nash, Anthony]

Dear all,


I understand that this is quite a basic question, but I think I need a
fresh set of eyes to figure out what¹s going on with gmx distance -select
syntax.

In my index file, I have 23 groups, two of which are ³ZINC² and ³CARBONYL²

Based on numerous mail archive suggestions and the gromacs website change
in tools page, I¹m using:

$CHEM/rungromacs505 gmx distance -n index.ndx -f umb_1_all.trr -s
umb_1_umb.tpr -oh umb_1_distance_hist.xvg -oall umb_1_distance.xvg -select
'com of group "ZINC" plus com of group "CARBONYL"'


Which kicks up the error:

Invalid command-line options
  In command-line option -select
In keyword 'com'
  required parameter 'of' not specified
invalid selection 'com'
  In command-line option -select
syntax error
invalid selection 'group'
  In command-line option -select
syntax error
invalid selection 'plus'
  In command-line option -select
In keyword 'com'
  required parameter 'of' not specified
invalid selection 'com'
  In command-line option -select
syntax error
invalid selection 'group'

I¹d appreciate any suggestions. I¹ve also tried to do *something* via the
interactive command line (when -select is not supplied), but what it¹s
asking for is a bit cryptic.

Thanks
Anthony  


-- 
Gromacs Users mailing list

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

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

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

Re: [gmx-users] Membrane protein insertion

[Nash, Anthony]

Hi,


I used this method recently and I was experiencing the same errors. As
Mark suggested, makes sure your protein survives an energy minimisation.
My error was a result of poor preparation of the .pdb file before running
pdb2gmx. My .itp file contained a long bond between the C and N termini of
a dimer. 

IF you have any doubts, try running InflateGRO2 using a very small part of
your protein (1 chain) and see if it works.

Anthony 

On 24/08/2015 08:59, Mark Abraham mark.j.abra...@gmail.com wrote:

Hi,

Something isn't stable. Check that your membrane protein survives a vacuum
EM (at least). And check your parameter settings for inflategro.

Mark

On Mon, Aug 24, 2015 at 9:53 AM Yasser Almeida Hernández 
yasser.almeida.hernan...@chemie.uni-hamburg.de wrote:

 Hi all,

 I want to insert a membrane protein into a model bilayer and I've tried
 several methods. The more straightforward seems to be Memembed, which
 gives an orientation to the membrane. On the other hand this method only
 outputs the membrane as dummy balls. How to translate this orientation
 to a model bilayer (DOPC, POPC)?

 On the other hand I've tried InflateGRO2 for this purpose and I got this
 error (following the tutorial in
 https://code.google.com/p/inflategro2/wiki/TutorialTolC ):

 Back Off! I just backed up ../1-topology/4hzuS_popc.top to
 ../1-topology/#4hzuS_popc.top.3#
 Will use a deflating factor of 0.925539817285086 to shrink the lipids
 back in 20 steps
 Deflating
 grompp -f deflate.mdp -c inflated.gro -p ../1-topology/4hzuS_popc.top -n
 inflated.ndx -o tmp.tpr -maxwarn 1
 mdrun -s tmp.tpr -v -deffnm tmp_out -c shrink.00.gro
-
 ERROR: Cannot open GRO file shrink.00.gro: No such file or directory


 The end of the log.out file is this:

 Fatal error:
 Too many LINCS warnings (1184)
 If you know what you are doing you can adjust the lincs warning
 threshold in your mdp file or set the environment variable
 GMX_MAXCONSTRWARN to -1, but normally it is better to fix the problem

 Any thoughts?

 Thanks in advance

 --
 Yasser Almeida Hernández
 PhD student
 Institute of Biochemistry and Molecular Biology
 Department of Chemistry
 University of Hamburg
 Martin-Luther-King-Platz 6
 20146 Hamburg
 Germany
 +49 40 42838 2845
 yasser.almeida.hernan...@chemie.uni-hamburg.de
 office: Grindelallee 117, room 250c

 --
 Gromacs Users mailing list

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

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

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

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

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

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

-- 
Gromacs Users mailing list

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

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

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

[gmx-users] Force field parameterisation

[Nash, Anthony]

Hi all,

I am trying to parameterise a new compound in Gromacs (theoretical
compound - my experimental collaborators are still trying to purify and
mass spec it) using the Amber 99sb force field. After asking about, I now
know that in Amber I would take my QM structure run antechamber (or R.E.D
tools) for partial charges and then parmchk2 for the bonded parameters. I
would then need to do *something* to get it in an .itp format.

I¹ve come a long way without need to use the Amber toolkit and I can¹t say
I want to start now (although I have just downloaded it). Can this be done
using anything Gromacs has to offer? Alternatively, does anyone know how
go from an Amber .lib file to an .itp file?

Many thanks
Anthony 


-- 
Gromacs Users mailing list

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

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

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

[gmx-users] High load imbalance: 31.8%

[Nash, Anthony]

Hi all,

I appear to have a very high load imbalance on some of my runs. Values
starting from approx. 7% up to 31.8% with reported vol min/aver of around
0.6 (I haven¹t found one under half yet).

When I look through the .log file at the start of the run I see:

Initializing Domain Decomposition on 8 ranks
Dynamic load balancing: auto
Will sort the charge groups at every domain (re)decomposition
Initial maximum inter charge-group distances:
two-body bonded interactions: 0.514 nm, LJ-14, atoms 3116 3123
  multi-body bonded interactions: 0.429 nm, Proper Dih., atoms 3116 3123
Minimum cell size due to bonded interactions: 0.472 nm
Maximum distance for 5 constraints, at 120 deg. angles, all-trans: 0.862 nm
Estimated maximum distance required for P-LINCS: 0.862 nm
This distance will limit the DD cell size, you can override this with -rcon
Using 0 separate PME ranks, as there are too few total
 ranks for efficient splitting
Scaling the initial minimum size with 1/0.8 (option -dds) = 1.25
Optimizing the DD grid for 8 cells with a minimum initial size of 1.077 nm
The maximum allowed number of cells is: X 12 Y 12 Z 12
Domain decomposition grid 4 x 2 x 1, separate PME ranks 0
PME domain decomposition: 4 x 2 x 1
Domain decomposition rank 0, coordinates 0 0 0
Using 8 MPI processes
Using 1 OpenMP thread per MPI process




Having a quick look through the documentation and I see that I should
consider implementing the verlet cut-off (which I am) and adjust the
number of PME nodes/cut-off and PME grid spacing. Would this simply be a
case of throwing more cores at the simulation or must I play around with
P-LINCS parameters?

Thanks
Anthony

-- 
Gromacs Users mailing list

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

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

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

Re: [gmx-users] High load imbalance: 31.8%

[Nash, Anthony]

Hi Mark,

Many thanks for looking into this.

One of the log files (the job hasn’t finished running) is here:
https://www.dropbox.com/s/zwrro54yni2uxtn/umb_3_umb.log?dl=0

The system is a soluble collagenase in water with a collagen substrate and
two zinc co-factors. There are 287562 atoms in the system.

Please let me know if you need to know anything else. Thanks!

Anthony





On 20/08/2015 11:39, Mark Abraham mark.j.abra...@gmail.com wrote:

Hi,

In cases like this, it's good to describe what's in your simulation, and
share the full .log file on a file-sharing service, so we can see both the
things mdrun reports early and late.

Mark

On Thu, Aug 20, 2015 at 8:22 AM Nash, Anthony a.n...@ucl.ac.uk wrote:

 Hi all,

 I appear to have a very high load imbalance on some of my runs. Values
 starting from approx. 7% up to 31.8% with reported vol min/aver of
around
 0.6 (I haven¹t found one under half yet).

 When I look through the .log file at the start of the run I see:

 Initializing Domain Decomposition on 8 ranks
 Dynamic load balancing: auto
 Will sort the charge groups at every domain (re)decomposition
 Initial maximum inter charge-group distances:
 two-body bonded interactions: 0.514 nm, LJ-14, atoms 3116 3123
   multi-body bonded interactions: 0.429 nm, Proper Dih., atoms 3116 3123
 Minimum cell size due to bonded interactions: 0.472 nm
 Maximum distance for 5 constraints, at 120 deg. angles, all-trans:
0.862 nm
 Estimated maximum distance required for P-LINCS: 0.862 nm
 This distance will limit the DD cell size, you can override this with
-rcon
 Using 0 separate PME ranks, as there are too few total
  ranks for efficient splitting
 Scaling the initial minimum size with 1/0.8 (option -dds) = 1.25
 Optimizing the DD grid for 8 cells with a minimum initial size of 1.077
nm
 The maximum allowed number of cells is: X 12 Y 12 Z 12
 Domain decomposition grid 4 x 2 x 1, separate PME ranks 0
 PME domain decomposition: 4 x 2 x 1
 Domain decomposition rank 0, coordinates 0 0 0
 Using 8 MPI processes
 Using 1 OpenMP thread per MPI process




 Having a quick look through the documentation and I see that I should
 consider implementing the verlet cut-off (which I am) and adjust the
 number of PME nodes/cut-off and PME grid spacing. Would this simply be a
 case of throwing more cores at the simulation or must I play around with
 P-LINCS parameters?

 Thanks
 Anthony

 --
 Gromacs Users mailing list

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

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

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

-- 
Gromacs Users mailing list

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

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

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

-- 
Gromacs Users mailing list

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

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

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

Re: [gmx-users] High load imbalance: 31.8%

[Nash, Anthony]

Hi Szilárd

Thanks for all of that advice. I’m going to have to take a lot of this up
with the Cluster Service Staff. This is a new cluster service I won a
grant for, thus not my usual platform which would typically yield an
imbalance of somewhere around 0.8% to 2%.

Thanks again
Anthony



On 20/08/2015 16:52, Szilárd Páll pall.szil...@gmail.com wrote:

Hi,

You're not pinning threads and it seems that you're running on a large SMP
machine! Assuming that the 512 threads reported (line 91) is correct
that's
a 32 socket SMP machine, perhaps an SGI UV? In any case Xeon E5-4xxx is
typically deployed in 4-8 socket installations, so your 8 threads will be
floating around on a number of CPUs which ruins your performance - and
likely contributes to the varying and large load imbalance.

My advice:
- don't ignore notes/warnings issued by mdrun (line 366, should be on the
standard out too), we put quite some though into spamming users only when
relevant :)
- pin mdrun and/or its threads either with -pin on (and -pinoffset if
needed) or with whatever tools your admins provide/recommend

[Extras: consider using FFTW even with the Intel compilers it's often
faster for our small FFTs than MKL; and GNU iso Intel compiler is often
faster too.]

Fixing the above issues should not only reduce imbalance but most likely
also allow you to gain quite some simulation performance! Let us know if
it
worked.

Cheers,

--
Szilárd

On Thu, Aug 20, 2015 at 5:08 PM, Nash, Anthony a.n...@ucl.ac.uk wrote:

 Hi Mark,

 Many thanks for looking into this.

 One of the log files (the job hasn’t finished running) is here:
 https://www.dropbox.com/s/zwrro54yni2uxtn/umb_3_umb.log?dl=0

 The system is a soluble collagenase in water with a collagen substrate
and
 two zinc co-factors. There are 287562 atoms in the system.

 Please let me know if you need to know anything else. Thanks!

 Anthony





 On 20/08/2015 11:39, Mark Abraham mark.j.abra...@gmail.com wrote:

 Hi,
 
 In cases like this, it's good to describe what's in your simulation,
and
 share the full .log file on a file-sharing service, so we can see both
the
 things mdrun reports early and late.
 
 Mark
 
 On Thu, Aug 20, 2015 at 8:22 AM Nash, Anthony a.n...@ucl.ac.uk wrote:
 
  Hi all,
 
  I appear to have a very high load imbalance on some of my runs.
Values
  starting from approx. 7% up to 31.8% with reported vol min/aver of
 around
  0.6 (I haven¹t found one under half yet).
 
  When I look through the .log file at the start of the run I see:
 
  Initializing Domain Decomposition on 8 ranks
  Dynamic load balancing: auto
  Will sort the charge groups at every domain (re)decomposition
  Initial maximum inter charge-group distances:
  two-body bonded interactions: 0.514 nm, LJ-14, atoms 3116 3123
multi-body bonded interactions: 0.429 nm, Proper Dih., atoms 3116
3123
  Minimum cell size due to bonded interactions: 0.472 nm
  Maximum distance for 5 constraints, at 120 deg. angles, all-trans:
 0.862 nm
  Estimated maximum distance required for P-LINCS: 0.862 nm
  This distance will limit the DD cell size, you can override this with
 -rcon
  Using 0 separate PME ranks, as there are too few total
   ranks for efficient splitting
  Scaling the initial minimum size with 1/0.8 (option -dds) = 1.25
  Optimizing the DD grid for 8 cells with a minimum initial size of
1.077
 nm
  The maximum allowed number of cells is: X 12 Y 12 Z 12
  Domain decomposition grid 4 x 2 x 1, separate PME ranks 0
  PME domain decomposition: 4 x 2 x 1
  Domain decomposition rank 0, coordinates 0 0 0
  Using 8 MPI processes
  Using 1 OpenMP thread per MPI process
 
 
 
 
  Having a quick look through the documentation and I see that I should
  consider implementing the verlet cut-off (which I am) and adjust the
  number of PME nodes/cut-off and PME grid spacing. Would this simply
be a
  case of throwing more cores at the simulation or must I play around
with
  P-LINCS parameters?
 
  Thanks
  Anthony
 
  --
  Gromacs Users mailing list
 
  * Please search the archive at
  http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
  posting!
 
  * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
 
  * For (un)subscribe requests visit
  https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
  send a mail to gmx-users-requ...@gromacs.org.
 
 --
 Gromacs Users mailing list
 
 * Please search the archive at
 http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
 posting!
 
 * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
 
 * For (un)subscribe requests visit
 https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
 send a mail to gmx-users-requ...@gromacs.org.

 --
 Gromacs Users mailing list

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

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

 * For (un)subscribe requests visit
 https://maillist.sys.kth.se

[gmx-users] distance issues with umbrella sampling

[Nash, Anthony]

As far as I¹ve understood, the absolute distance reported using g_dist
(note: alternative name in 5+) and the reported harmonic potential between
two groups using pull code in grompp, doesn¹t always match.

As such, I some times end up with neighbouring umbrella histograms
practically sitting on top of one another rather than nicely overlapping
at the tails. Would anyone with more experience than myself at umbrella
sampling clarify whether I should keep in both histograms, or remove one
of the hisograms. 

Also, if my reaction coordinate is an absolute distance and therefor I
have set pull_dim = Y Y Y (rather than, say, Z for water through an ion
channel or Y for TM protein dimerisation), is it appropriate to use the
pullf.xvg and the window.tpr files in g_wham? Or is this beyond the
capabilities of g_wham. You see, my reaction coordinate is still distance
but it was defined via a traversed trajectory over a potential energy
landscape using meta-dynamics.

Many thanks
Anthony  

Dr Anthony Nash
Department of Chemistry
University College London


-- 
Gromacs Users mailing list

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

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

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

Re: [gmx-users] Change to umbrella sampling pull code

[Nash, Anthony]

Thanks for clearing that up Justin, I’ve adjusted accordingly. All but one
error remains:


Pull group 1 'ZINC' has 1 atoms
Pull group 2 'CARBONYL' has 1 atoms

---
Program grompp, VERSION 5.0.5
Source code file: 
/home/columbus/chem_software/GROMACS/2015/gromacs-5.0.5/src/gromacs/gmxprep
rocess/readpull.c, line: 377


Fatal error:
Identical pull group indices in pull-coord1-groups
For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors
---


The entries in my index.ndx are:

[ ZINC ]
5809
[ CARBONYL ]
6481


For the respective .mdp entries (pull code):

; Pull code
pull= umbrella
pull_geometry   = distance
pull_dim= Y Y Y
pull_start  = yes
pull_ngroups= 2
pull_group1_name = ZINC
pull_group2_name = CARBONYL
pull_coord1_init = 0
pull_coord1_rate = 0
pull_coord1_k = 1000
pull_nstxout= 1000  ; every 2 ps
pull_nstfout= 1000  ; every 2 ps




On 12/08/2015 21:26, Justin Lemkul jalem...@vt.edu wrote:



On 8/12/15 3:20 PM, Nash, Anthony wrote:
 Hi Mark,

 Thanks for the reply.

 Comparing the documentation between manual-4.5.6 (what I had previously
 been using) and manual-5.0.4:

 4.5.6 has reference to pull_group0 as reference group, and in 5.0.4
there
 is no explicit mention (that I can see) to a reference group but there
is
 one for the pull group (pull-group1-name).


Right, because multiple reaction coordinates can now be simultaneously be
restrained.  There is no requirement for a single reference group.  Any
combination of groups can be used.  Offhand, your settings are right
except 
pull-ngroups is 2, as there are two groups.

-Justin

 Thanks
 Anthony

 Dr Anthony Nash
 Department of Chemistry
 University College London





 On 12/08/2015 19:35, Mark Abraham mark.j.abra...@gmail.com wrote:

 Hi,
 On Wed, Aug 12, 2015 at 6:20 PM Nash, Anthony a.n...@ucl.ac.uk wrote:

 Dear all,

 This is the first time I¹ve ran pull code (for umbrella sampling)
since
 the change from Gromacs 4 to gromancs 5. I¹ve noticed some difference
in
 the .mdp key-value parameters. Could I have a sanity check on the
values
 below.


 I have no real idea, but you could start by comparing the respective
 documentation. If something's not clear, that's something to fix ;-)

 Also, given that I want a harmonic potential between the two
 groups, does it matter which index group is assigned to the respective
 Œpull-groupN-name¹ key?


 I doubt it, but if you want a definitive answer, make yourself a toy
test
 case and see.

 Mark


 Note: my reaction coordinate was defined initially following the
 trajectory between two minima during a meta-dynamics MD simulation.
That
 too was distance based.


 ;Pull code
 pull= umbrella
 pull-geometry   = distance
 pull-dim= Y Y Y
 pull-start  = yes
 pull-ngroups= 1
 pull-group1-name = ZINC
 pull-group2-name = CARBONYL
 pull-coord1-init = 0
 pull-coord1-rate = 0
 pull-coord1-k = 1000
 pull-nstxout= 1000  ; every 2 ps
 pull-nstfout= 1000  ; every 2 ps


 Many thanks
 Anthony

 --
 Gromacs Users mailing list

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

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

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

 --
 Gromacs Users mailing list

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

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

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


-- 
==

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

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

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

==
-- 
Gromacs Users mailing list

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

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

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

-- 
Gromacs Users mailing list

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

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

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman

[gmx-users] Change to umbrella sampling pull code

[Nash, Anthony]

Dear all,

This is the first time I¹ve ran pull code (for umbrella sampling) since
the change from Gromacs 4 to gromancs 5. I¹ve noticed some difference in
the .mdp key-value parameters. Could I have a sanity check on the values
below. Also, given that I want a harmonic potential between the two
groups, does it matter which index group is assigned to the respective
Œpull-groupN-name¹ key?

Note: my reaction coordinate was defined initially following the
trajectory between two minima during a meta-dynamics MD simulation. That
too was distance based.


;Pull code
pull= umbrella
pull-geometry   = distance
pull-dim= Y Y Y
pull-start  = yes
pull-ngroups= 1
pull-group1-name = ZINC
pull-group2-name = CARBONYL
pull-coord1-init = 0
pull-coord1-rate = 0
pull-coord1-k = 1000
pull-nstxout= 1000  ; every 2 ps
pull-nstfout= 1000  ; every 2 ps


Many thanks
Anthony

-- 
Gromacs Users mailing list

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

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

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

Re: [gmx-users] Change to umbrella sampling pull code

[Nash, Anthony]

Hi Mark,

Thanks for the reply.

Comparing the documentation between manual-4.5.6 (what I had previously
been using) and manual-5.0.4:

4.5.6 has reference to pull_group0 as reference group, and in 5.0.4 there
is no explicit mention (that I can see) to a reference group but there is
one for the pull group (pull-group1-name).

Thanks
Anthony

Dr Anthony Nash
Department of Chemistry
University College London





On 12/08/2015 19:35, Mark Abraham mark.j.abra...@gmail.com wrote:

Hi,
On Wed, Aug 12, 2015 at 6:20 PM Nash, Anthony a.n...@ucl.ac.uk wrote:

 Dear all,

 This is the first time I¹ve ran pull code (for umbrella sampling) since
 the change from Gromacs 4 to gromancs 5. I¹ve noticed some difference in
 the .mdp key-value parameters. Could I have a sanity check on the values
 below.


I have no real idea, but you could start by comparing the respective
documentation. If something's not clear, that's something to fix ;-)

Also, given that I want a harmonic potential between the two
 groups, does it matter which index group is assigned to the respective
 Œpull-groupN-name¹ key?


I doubt it, but if you want a definitive answer, make yourself a toy test
case and see.

Mark


 Note: my reaction coordinate was defined initially following the
 trajectory between two minima during a meta-dynamics MD simulation. That
 too was distance based.


 ;Pull code
 pull= umbrella
 pull-geometry   = distance
 pull-dim= Y Y Y
 pull-start  = yes
 pull-ngroups= 1
 pull-group1-name = ZINC
 pull-group2-name = CARBONYL
 pull-coord1-init = 0
 pull-coord1-rate = 0
 pull-coord1-k = 1000
 pull-nstxout= 1000  ; every 2 ps
 pull-nstfout= 1000  ; every 2 ps


 Many thanks
 Anthony

 --
 Gromacs Users mailing list

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

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

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

-- 
Gromacs Users mailing list

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

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

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

-- 
Gromacs Users mailing list

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

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

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

[gmx-users] Umbrella sampling sanity check

[Nash, Anthony]

Hi all,

I would appreciate a little sanity check for umbrella sampling pull code
parameters.

My reaction coordinate is defined as the distance between a globular
soluble enzyme and a substrate in solution (water). I have already
captured the individual windows using the trajectory generated from an
earlier meta-dynamics simulation. I would now like to perform umbrella
sampling, a technique which I feel is more accountable for the free
energies it generates once a reaction coordinate is known.

My two questions:
1) My earlier meta-dynamics simulation used the distance between the
enzyme active site zinc ion to a particular substrate atom. The substrate
atom (a backbone carbonyl oxygen) would in reality coordinate with the
zinc prior to proteolysis. In pull code, can pull_group0 and pull_group1
simply be the same two atoms (Zn2+ and O)? As the distance harmonic
potential works on COM, this would just be a point of an atoms.
 
2) Up until now I¹ve always worked on bilayer, so forgive the simplicity.
Given that my substrate could diffuse in any direction, should I set
pull_dim=Y Y Y?

Many thanks
Anthony

Dr Anthony Nash
Department of Chemistry
University College London





-- 
Gromacs Users mailing list

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

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

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

Re: [gmx-users] Output frames within a defined distance range

[Nash, Anthony]

Hi Justin,

Cheers, that was perfect.

In case anyone comes across this in months from now, I simply performed a
gmx_d distance command which generates a .xvg of [timestamp] [distance].
This is then used as the input -drop for gmx_d trjconv. Combining with
-dropunder [mindistance] and -dropover [maxdistance], you end up with a
.trr/.xtc/.gro file with only the frames that meets your distance
criteria. 

Thanks
Anthony 

Dr Anthony Nash
Department of Chemistry
University College London





On 06/08/2015 13:03, Justin Lemkul jalem...@vt.edu wrote:



On 8/6/15 3:56 AM, Nash, Anthony wrote:
 Hi all,

 I¹m hoping to avoid putting together a simple script. Is there an option
 in Gromacs 5+ to only output frames from a trajectory if a certain
 criteria is meet? In this case, I only want the frames from a trajectory
 if the carbon-alpha of a particular GLY residue is within 0.8 - 1 nm of
a
 zinc ion.


I think this is what gmx trjconv -drop -dropover/-dropunder does.  I have
never 
tried it, though.

-Justin

-- 
==

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

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

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

==
-- 
Gromacs Users mailing list

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

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

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

-- 
Gromacs Users mailing list

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

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

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

[gmx-users] Output frames within a defined distance range

[Nash, Anthony]

Hi all,

I¹m hoping to avoid putting together a simple script. Is there an option
in Gromacs 5+ to only output frames from a trajectory if a certain
criteria is meet? In this case, I only want the frames from a trajectory
if the carbon-alpha of a particular GLY residue is within 0.8 - 1 nm of a
zinc ion. 

Many thanks
Anthony 


-- 
Gromacs Users mailing list

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

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

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

Re: [gmx-users] Enforced rotation errors

[Nash, Anthony]

Dear Carsten,

Thanks for that suggestion. Seems like that solved that particular
problem. Unfortunately though, my trajectory is not what I expected.

I have been able to rotate the cylindrical ligand about it¹s principle
axis in an earlier Œprototype¹ of my enzyme-ligand complex using iso-pf
enforced rotation. But after replacing the ligand with the Œreal¹ thing
and running a regular MD simulation for 200ns, I¹ve failed to recreate the
desired rotation using enforced rotation. I original started with iso-pf
(as it worked in the earlier case), but this resulted in lots of LINCS
errors. I then tried flex (as per my original post) and now flex-t.

I begin by calculating the moment of inertia of my cylindrical enzyme.
Using the lal01psx script for VMD I pick the 3rd vector which I add to the
.mdp file as:

rot-vec0 = 0.691585601358247 -0.6995947474399045 0.17965674311989613

It looks as though the ligand is being rotated in a big arch, rather than
rotated around its length/long axis. I take it enforced rotation, applies
the potential about this vector? My inputs are:

; ENFORCED ROTATION
; Enforced rotation: No or Yes
rotation = Yes
; Output frequency for angle, torque and rotation potential energy for the
whole group
rot-nstrout  = 1
; Output frequency for per-slab data (angles, torques and slab centers)
rot-nstsout  = 10
; Number of rotation groups
rot-ngroups  = 1
; Rotation group name
rot-group0   = Collagen_CA
; Rotation potential. Can be iso, iso-pf, pm, pm-pf, rm, rm-pf, rm2,
rm2-pf, flex, flex-t, flex2, flex2-t
rot-type0= flex-t ;iso-pf ;using a pivot free i.e., a
detached!:
; Use mass-weighting of the rotation group positions
rot-massw0   = yes
; Rotation vector, will get normalized
rot-vec0 = 0.691585601358247 -0.6995947474399045 0.17965674311989613


; Pivot point for the potentials iso, pm, rm, and rm2 [nm]
;rot-pivot0   = 4.31852e+00  1.73201e+00  1.89800e+00


; Rotation rate [degree/ps] and force constant [kJ/(mol*nm^2)]
rot-rate0= 0.021  ;
rot-k0   = 600.0 ; - change 100 through to 600
; Slab distance for flexible axis rotation [nm]
rot-slab-dist0   = 1
; Minimum value of Gaussian function for the force to be evaluated (for
flex* potentials)
rot-min-gauss0   = 0.001
; Value of additive constant epsilon' [nm^2] for rm2* and flex2* potentials
rot-eps0 = 0.0001
; Fitting method to determine angle of rotation group (rmsd, norm, or
potential)
rot-fit-method0  = norm
; For fit type 'potential', nr. of angles around the reference for which
the pot. is evaluated
rot-potfit-nsteps0   = 21
; For fit type 'potential', distance in degrees between two consecutive
angles
rot-potfit-step0 = 0.25

Any suggestions are gratefully appreciated.

Thanks
Anthony

Dr Anthony Nash
Department of Chemistry
University College London





On 20/07/2015 15:53, Kutzner, Carsten ckut...@gwdg.de wrote:

Dear Anthony,

the problem you are experiencing with the Œflex¹ rotation potential
could be related to the rotation group moving too far along the direction
of the rotation vector. As for V_flex, the slabs are fixed in space,
the rotation group may after some time enter a region where no reference
slab centers are defined, triggering the error that you see.

In that case, you can use the Œflex-t¹ variant of the potential. Here,
the midplane of slab n=0 is attached to the center of mass of the
rotation group, so that the slabs move with the rotation group.
See equations 6.46 and 6.47 in the GROMACS 5.0 PDF manual.

Carsten


 On 20 Jul 2015, at 16:20, Nash, Anthony a.n...@ucl.ac.uk wrote:
 
 Dear All,
 
 I hope you can help. I am using 'flex' enforced rotation to rotate a
cylindrical protein along the surface of a globular protein.
Unfortunately my system is experiencing what I can only assume is an IO
problem:
 
 DD  step 94 load imb.: force  2.9%  pme mesh/force 0.677
 
   Step   Time Lambda
 95 1900.00.0
 
   Energies (kJ/mol)
  AngleProper Dih.  Improper Dih.  LJ-14
Coulomb-14
1.57412e+041.99606e+049.25263e+027.98071e+03
8.61345e+04
LJ (SR)  Disper. corr.   Coulomb (SR)   Coul. recip. Position
Rest.
9.59944e+05   -8.46602e+04   -6.08542e+066.95493e+04
1.22242e+03
   COM Pull En.  PotentialKinetic En.   Total Energy
Temperature
3.91316e+02   -5.00823e+068.56761e+05   -4.15147e+06
3.10369e+02
 Pres. DC (bar) Pressure (bar)   Constr. rmsd
   -4.21942e+023.27704e+002.75540e-05
 
 
 ---
 Program mdrun_mpi, VERSION 5.0.5
 Source code file:
/home/columbus/chem_software/GROMACS/2015/gromacs-5.0.5/src/gromacs/pulli
ng/pull_rotation.c, line: 2502
 
 Fatal error:
 Enforced rotation: No reference data for first slab (n=-13), unable to
proceed.
 For more

[gmx-users] Enforced rotation errors

[Nash, Anthony]

Dear All,

I hope you can help. I am using 'flex' enforced rotation to rotate a 
cylindrical protein along the surface of a globular protein. Unfortunately my 
system is experiencing what I can only assume is an IO problem:

DD  step 94 load imb.: force  2.9%  pme mesh/force 0.677

   Step   Time Lambda
 95 1900.00.0

   Energies (kJ/mol)
  AngleProper Dih.  Improper Dih.  LJ-14 Coulomb-14
1.57412e+041.99606e+049.25263e+027.98071e+038.61345e+04
LJ (SR)  Disper. corr.   Coulomb (SR)   Coul. recip. Position Rest.
9.59944e+05   -8.46602e+04   -6.08542e+066.95493e+041.22242e+03
   COM Pull En.  PotentialKinetic En.   Total EnergyTemperature
3.91316e+02   -5.00823e+068.56761e+05   -4.15147e+063.10369e+02
 Pres. DC (bar) Pressure (bar)   Constr. rmsd
   -4.21942e+023.27704e+002.75540e-05


---
Program mdrun_mpi, VERSION 5.0.5
Source code file: 
/home/columbus/chem_software/GROMACS/2015/gromacs-5.0.5/src/gromacs/pulling/pull_rotation.c,
 line: 2502

Fatal error:
Enforced rotation: No reference data for first slab (n=-13), unable to proceed.
For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors
---


I am using gromacs 5.0.5. The .mdp settings are:

; ENFORCED ROTATION
; Enforced rotation: No or Yes
rotation = Yes
; Output frequency for angle, torque and rotation potential energy for the 
whole group
rot-nstrout  = 1
; Output frequency for per-slab data (angles, torques and slab centers)
rot-nstsout  = 10
; Number of rotation groups
rot-ngroups  = 1
; Rotation group name
rot-group0   = Collagen_CA
; Rotation potential. Can be iso, iso-pf, pm, pm-pf, rm, rm-pf, rm2, rm2-pf, 
flex, flex-t, flex2, flex2-t
rot-type0= flex
; Use mass-weighting of the rotation group positions
rot-massw0   = yes
; Rotation vector, will get normalized
rot-vec0 = 0.691585601358247 -0.6995947474399045 0.17965674311989613

; Pivot point for the potentials iso, pm, rm, and rm2 [nm]
;rot-pivot0   = 4.31852e+00  1.73201e+00  1.89800e+00

; Rotation rate [degree/ps] and force constant [kJ/(mol*nm^2)]
rot-rate0= 0.021
rot-k0   = 600.0 ; testing these
; Slab distance for flexible axis rotation [nm]
rot-slab-dist0   = 1
; Minimum value of Gaussian function for the force to be evaluated (for flex* 
potentials)
rot-min-gauss0   = 0.001
; Value of additive constant epsilon' [nm^2] for rm2* and flex2* potentials
rot-eps0 = 0.0001
; Fitting method to determine angle of rotation group (rmsd, norm, or potential)
rot-fit-method0  = norm
; For fit type 'potential', nr. of angles around the reference for which the 
pot. is evaluated
rot-potfit-nsteps0   = 21
; For fit type 'potential', distance in degrees between two consecutive angles
rot-potfit-step0 = 0.25

The following files are being generated:

flex_1.trr
flex_1.edr
flex_1.log
flex_1.xvg
#flex_1.log.1#
#flex_1.log.2#
#flex_1.log.3#
flex_1_out   --output from the cluster

Any help is greatly appreciated.
Kind regards
Anthony

Dr Anthony Nash
Department of Chemistry
University College London

-- 
Gromacs Users mailing list

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

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

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

[gmx-users] Gromacs with InflateGro

[Nash, Anthony]

Hi all,

I'm really struggling to get InflateGro to work in a dimer + lipid (no water or 
ions) system. The fault seems to happen when the regular expression to break 
the .gro entry reads in an entry from C  to C1. I managed to generalise 
the regex further and now using substr to explicitly pull out the entries I 
want. Unfortunately I'm a little confused as to the purpose of a variable 
called 'now'.  Has anyone got a bullet proof copy of the code working with 
+ atoms as I've had enough going through this file.

I also tried InflateGro2, which passes the first two iterations but then dies 
after the second. It can't find tmp_out.gro, as the previous entry hadn't 
successfully passed grompp due to a difference in atom count between .gro and 
.top. 


Thanks
Anthony 
-- 
Gromacs Users mailing list

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

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

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

[gmx-users] grompp_d fatal error in Amber FF

[Nash, Anthony]

Hi all,

Gromacs ver 5.0.4 - no forcefield modification (although I had, but changed 
back for testing purposes). Installed on a brand new machine: this is my 
laptop's first grompp. 

I was trying to inflate and deflate a bilayer using InflateGro but it through 
up a fatal error. I removed everything to do with that code and ran a regular 
grompp with a very sparse .mdp file.


I'm getting the following error:

Program grompp_d, VERSION 5.0.4
Source code file: 
/Users/acnash/Gromacs/gromacs-5.0.4/src/gromacs/gmxpreprocess/topio.c, line: 633

Fatal error:
Unknown error - File amber99sb-ildn.ff, line 0
Last line read:
' '

The matching code for this is:
do
 626 {
 627 status = cpp_read_line(handle, STRLEN, line);
 628 done   = (status == eCPP_EOF);
 629 if (!done)
 630 {
 631 if (status != eCPP_OK)
 632 {
 633 gmx_fatal(FARGS, cpp_error(handle, status));
 634 }

So obviously it hasn't read something in.

Anthony thoughts.

Thanks
Anthony



-- 
Gromacs Users mailing list

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

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

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

Re: [gmx-users] grompp_d fatal error in Amber FF

[Nash, Anthony]

Hi All,

My apologies, false alarm. 

The issue was with specifying the forcefield FILE (NOT the directory) in my 
topology file. 

Thanks
Anthony

Dr Anthony Nash
Department of Chemistry
University College London

From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se 
[gromacs.org_gmx-users-boun...@maillist.sys.kth.se] on behalf of Nash, Anthony 
[a.n...@ucl.ac.uk]
Sent: 07 March 2015 08:05
To: gmx-us...@gromacs.org
Subject: [gmx-users]  grompp_d fatal error in Amber FF

Hi all,

Gromacs ver 5.0.4 - no forcefield modification (although I had, but changed 
back for testing purposes). Installed on a brand new machine: this is my 
laptop's first grompp.

I was trying to inflate and deflate a bilayer using InflateGro but it through 
up a fatal error. I removed everything to do with that code and ran a regular 
grompp with a very sparse .mdp file.


I'm getting the following error:

Program grompp_d, VERSION 5.0.4
Source code file: 
/Users/acnash/Gromacs/gromacs-5.0.4/src/gromacs/gmxpreprocess/topio.c, line: 633

Fatal error:
Unknown error - File amber99sb-ildn.ff, line 0
Last line read:
' '

The matching code for this is:
do
 626 {
 627 status = cpp_read_line(handle, STRLEN, line);
 628 done   = (status == eCPP_EOF);
 629 if (!done)
 630 {
 631 if (status != eCPP_OK)
 632 {
 633 gmx_fatal(FARGS, cpp_error(handle, status));
 634 }

So obviously it hasn't read something in.

Anthony thoughts.

Thanks
Anthony



--
Gromacs Users mailing list

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

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

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

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

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

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

[gmx-users] g_sas values of zero and Warning

[Nash, Anthony]

HI All,

I'm trying to get a measure of the solvent accessible surface area of a 
protein's catalytic site. It is unknown precisely how the substrate actually 
fits in the site given that in the crystal structure the site it too enclosed 
for the bulky substrate. I am using a progressive measure of g_sas -probe to 
get an idea of how accessible the catalytic site is during the simulation.

So I've tried six different -probe sizes, 0.1 (nm) 0.8, default , 1.6, 2.4, 
3.2. Only in the first two instances do I actually get any read out. It appears 
as though I am unable to access the active site residues as the probe radius is 
near to a vdw radius. I could asume that the catalytic site is very enclosed, 
which could physiologically make sense. Or I have done something wrong, after 
all I am getting the warning:


WARNING: could not find a Van der Waals radius for 2 atoms
205 out of 205 atoms were classified as hydrophobic

Generated values per probe radius below:


-probe 0.1 (nm)
 0 8.30466   0 8.30466
   100 8.04708   0 8.04708
   200 7.82871   0 7.82871
.

DEFAULT: -probe   real   0.14Radius of the solvent probe (nm)
 0 4.30913   0 4.30913
   100 4.24779   0 4.24779
   200 3.81957   0 3.81957

-probe 0.8 (nm)
 0   0   0   0
   100   0   0   0
   200   0   0   0
..

-probe 1.6 (nm)
 0   0   0   0
   100   0   0   0
   200   0   0   0
..

-probe 2.4
 0   0   0   0
   100   0   0   0
   200   0   0   0


-probe 3.2
 0   0   0   0
   100   0   0   0
   200   0   0   0

Many thanks
Anthony
-- 
Gromacs Users mailing list

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

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

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

Re: [gmx-users] Changing number of processors after a job restart

[Nash, Anthony]

Hi Justin,

I thought I would just give it a try at the risk of then being dropped back 
down on the queue. Turns out I end up with warnings rather than system-halting 
errors:

-
  #nodes mismatch,
current program: 32
checkpoint file: 64

  #PME-nodes mismatch,
current program: -1
checkpoint file: 16

Gromacs binary or parallel settings not identical to previous run.
Continuation is exact, but is not guaranteed to be binary identical.

--

I am now weighing up the pros and cons of this simulation. 

Thanks
Anthony

Dr Anthony Nash
Department of Chemistry
University College London

From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se 
[gromacs.org_gmx-users-boun...@maillist.sys.kth.se] on behalf of Justin Lemkul 
[jalem...@vt.edu]
Sent: 09 January 2015 12:40
To: gmx-us...@gromacs.org
Subject: Re: [gmx-users] Changing number of processors after a job restart

On 1/9/15 2:00 AM, Nash, Anthony wrote:
 Hi all,

 This is probably quite a fundamental bit of knowledge I am missing (and 
 struggling to find). In an effort to just get a system running rather than 
 waiting on a queue I am considering taking my job which has already ran for 
 48 hours and reducing the requested number of nodes. I would use the usual 
 -cpi .cpt -noappend notation in the job script to resubmit.

 I have a feeling though, that all manor of parallel calculations were 
 preserved in the check point file and are loaded upon restart. Would my job 
 reload and recalculate all the relevant cell decomposition, etc., without 
 throwing up an error.


IIRC, you will get errors when reading from the checkpoint file.  The
alternative is to re-create the .tpr with grompp (using -t to pass the .cpt)
rather than the usual tpbconv/mdrun -cpi approach.

-Justin

--
==

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

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

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

==
--
Gromacs Users mailing list

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

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

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

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

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

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

Re: [gmx-users] Changing number of processors after a job restart

[Nash, Anthony]

Fantastic! I figured it would all fall under the dominion of floating point 
arithmetic reproducibility.

Thanks
Anthony 

Dr Anthony Nash
Department of Chemistry
University College London

From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se 
[gromacs.org_gmx-users-boun...@maillist.sys.kth.se] on behalf of Justin Lemkul 
[jalem...@vt.edu]
Sent: 09 January 2015 12:49
To: gmx-us...@gromacs.org
Subject: Re: [gmx-users] Changing number of processors after a job restart

On 1/9/15 7:46 AM, Nash, Anthony wrote:
 Hi Justin,

 I thought I would just give it a try at the risk of then being dropped back 
 down on the queue. Turns out I end up with warnings rather than 
 system-halting errors:

 -
#nodes mismatch,
  current program: 32
  checkpoint file: 64

#PME-nodes mismatch,
  current program: -1
  checkpoint file: 16

 Gromacs binary or parallel settings not identical to previous run.
 Continuation is exact, but is not guaranteed to be binary identical.

 --

 I am now weighing up the pros and cons of this simulation.


Any possible differences are unlikely to be more significant than the normal
chaos of MD, I would think.  See, for instance:

http://www.gromacs.org/Documentation/How-tos/Extending_Simulations#Exact_vs_binary_identical_continuation

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

-Justin

--
==

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

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

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

==
--
Gromacs Users mailing list

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

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

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

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

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

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

[gmx-users] Changing number of processors after a job restart

[Nash, Anthony]

Hi all,

This is probably quite a fundamental bit of knowledge I am missing (and 
struggling to find). In an effort to just get a system running rather than 
waiting on a queue I am considering taking my job which has already ran for 48 
hours and reducing the requested number of nodes. I would use the usual -cpi 
.cpt -noappend notation in the job script to resubmit.

I have a feeling though, that all manor of parallel calculations were preserved 
in the check point file and are loaded upon restart. Would my job reload and 
recalculate all the relevant cell decomposition, etc., without throwing up an 
error.

Many thanks
Anthony 
-- 
Gromacs Users mailing list

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

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

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

Re: [gmx-users] Umbrella sampling along a dihedral angle

[Nash, Anthony]

Hi Justin,

Thanks for the reply, always appreciated. I like the way gromacs can handle 
modular structures by breaking up the topology into separate .itp file. 
Unfortunately this is the first time I've had to span two subunits with a 
specific restraint. I quickly knocked together a perl script to combine .itp 
topology files, and now it's all working fine. 

Thanks again. 
Anthony

Dr Anthony Nash
Department of Chemistry
University College London

From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se 
[gromacs.org_gmx-users-boun...@maillist.sys.kth.se] on behalf of Justin Lemkul 
[jalem...@vt.edu]
Sent: 28 December 2014 19:57
To: gmx-us...@gromacs.org
Subject: Re: [gmx-users] Umbrella sampling along a dihedral angle

On 12/27/14 1:50 PM, Nash, Anthony wrote:
 Hi Justin,

 I think I've shot myself in the foot and I'm trawling the gromacs manual with 
 little success at the moment. I'm trying to define my [ dihedral_restraint ] 
 over my ligand-enzyme complex. The topology of the ligand is in one .itp file 
 and the topology of the enzyme is in a separate .itp file. Therefore, in my 
 top level topology file (system.top) I have:

 #include ../enzyme.itp
 #include ../ligand.itp

 Typically if I want restraints to help equilibrate the system I would just 
 define a posre_NAME.itp with the respective force constants beneath each 
 #include directive. But given that each topology file starts with atom 1 I am 
 a bit unsure how I can define a set of restraints over two topology files.

 Can this be done? I'm really hoping I won't need to clump all topologies into 
 one big file.


Restraints are always numbered with respect to the [moleculetype] to which they
apply.  As an example:

http://www.gromacs.org/Documentation/Errors#Atom_index_n_in_position_restraints_out_of_bounds

If your (pseudo)dihedral spans the protein and the ligand, the only way to
proceed with a restraint is to have a merged [moleculetype] of protein and
ligand.  They simply cannot be in separate topologies in this case.

-Justin

 Thanks as ever,
 Anthony



 Dr Anthony Nash
 Department of Chemistry
 University College London
 
 From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se 
 [gromacs.org_gmx-users-boun...@maillist.sys.kth.se] on behalf of Justin 
 Lemkul [jalem...@vt.edu]
 Sent: 27 December 2014 14:49
 To: gmx-us...@gromacs.org
 Subject: Re: [gmx-users] Umbrella sampling along a dihedral angle

 On 12/26/14 6:11 PM, Nash, Anthony wrote:
 Dear Gromacs community,

 Using enforced rotation potentials I have generated a really smooth rotation
 of my ligand the catalytic binding domain of my protein. I then put together
 a simple perl script to calculate the dihedral angle of four particular atoms
 at each time step using g_angle. Taking a 6 degree interval (for a total of
 60 windows to span the 360 degree rotation) I modified the topology with a [
 dihedral_restraint ] with the respective dihedral angle for that particular
 configuation. I plan on then running each window for 10 ns to begin with
 (according to the PMF profile I can add more time). It is from this point
 where I am a little stuck.

 I have gone through a number of posts going back to 2011 with regards to
 using umbrella sampling along a dihedral angle reaction coordinate. The
 precise post was from Justin:

 You'll have to build various configurations that correspond to different
 dihedral angles (which form the sampling windows), then restrain them.

 The energy attributed to the restraints is then stored in the .edr file. From
 these energies, you should be able to construct the energy curve over the
 sampling windows. There are examples of this in the literature, so I suspect
 you should be able to find some demonstrations of how it's applied.

 I was wondering whether gromacs now support this procedure natively? I
 noticed g_wham has an option -cycl for dihedrals but I am still not sure how
 this fits in with the traditional gromacs means of generating pull force
 output as I have only been using distance based rather than angular based
 umbrella sampling.


 g_wham is still only really designed to process the output from the pull code,
 not any generalized restraint.  That would be a nice features, but AFAIK it is
 not in the pipeline.  You can de-bias the results like I suggested above, but
 it's probably more work than it's worth.  I would suggest just using Alan
 Grossfield's WHAM implementation.  There, you just need to provide the raw
 numbers (in this case, the dihedral values) and the force constants used.  It
 will then give you your PMF.  We have done this recently; it is very easy.

 -Justin

 --
 ==

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

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

Re: [gmx-users] Umbrella sampling along a dihedral angle

[Nash, Anthony]

Thanks for the help Justin. I've found the site: 

http://membrane.urmc.rochester.edu/content/wham

Cheers
Anthony 


Dr Anthony Nash
Department of Chemistry
University College London

From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se 
[gromacs.org_gmx-users-boun...@maillist.sys.kth.se] on behalf of Justin Lemkul 
[jalem...@vt.edu]
Sent: 27 December 2014 14:49
To: gmx-us...@gromacs.org
Subject: Re: [gmx-users] Umbrella sampling along a dihedral angle

On 12/26/14 6:11 PM, Nash, Anthony wrote:
 Dear Gromacs community,

 Using enforced rotation potentials I have generated a really smooth rotation
 of my ligand the catalytic binding domain of my protein. I then put together
 a simple perl script to calculate the dihedral angle of four particular atoms
 at each time step using g_angle. Taking a 6 degree interval (for a total of
 60 windows to span the 360 degree rotation) I modified the topology with a [
 dihedral_restraint ] with the respective dihedral angle for that particular
 configuation. I plan on then running each window for 10 ns to begin with
 (according to the PMF profile I can add more time). It is from this point
 where I am a little stuck.

 I have gone through a number of posts going back to 2011 with regards to
 using umbrella sampling along a dihedral angle reaction coordinate. The
 precise post was from Justin:

 You'll have to build various configurations that correspond to different
 dihedral angles (which form the sampling windows), then restrain them.

 The energy attributed to the restraints is then stored in the .edr file. From
 these energies, you should be able to construct the energy curve over the
 sampling windows. There are examples of this in the literature, so I suspect
 you should be able to find some demonstrations of how it's applied.

 I was wondering whether gromacs now support this procedure natively? I
 noticed g_wham has an option -cycl for dihedrals but I am still not sure how
 this fits in with the traditional gromacs means of generating pull force
 output as I have only been using distance based rather than angular based
 umbrella sampling.


g_wham is still only really designed to process the output from the pull code,
not any generalized restraint.  That would be a nice features, but AFAIK it is
not in the pipeline.  You can de-bias the results like I suggested above, but
it's probably more work than it's worth.  I would suggest just using Alan
Grossfield's WHAM implementation.  There, you just need to provide the raw
numbers (in this case, the dihedral values) and the force constants used.  It
will then give you your PMF.  We have done this recently; it is very easy.

-Justin

--
==

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

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

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

==
--
Gromacs Users mailing list

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

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

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

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

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

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

Re: [gmx-users] Umbrella sampling along a dihedral angle

[Nash, Anthony]

Hi Justin,

I think I've shot myself in the foot and I'm trawling the gromacs manual with 
little success at the moment. I'm trying to define my [ dihedral_restraint ] 
over my ligand-enzyme complex. The topology of the ligand is in one .itp file 
and the topology of the enzyme is in a separate .itp file. Therefore, in my top 
level topology file (system.top) I have:

#include ../enzyme.itp
#include ../ligand.itp 

Typically if I want restraints to help equilibrate the system I would just 
define a posre_NAME.itp with the respective force constants beneath each 
#include directive. But given that each topology file starts with atom 1 I am a 
bit unsure how I can define a set of restraints over two topology files. 

Can this be done? I'm really hoping I won't need to clump all topologies into 
one big file. 

Thanks as ever,
Anthony



Dr Anthony Nash
Department of Chemistry
University College London

From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se 
[gromacs.org_gmx-users-boun...@maillist.sys.kth.se] on behalf of Justin Lemkul 
[jalem...@vt.edu]
Sent: 27 December 2014 14:49
To: gmx-us...@gromacs.org
Subject: Re: [gmx-users] Umbrella sampling along a dihedral angle

On 12/26/14 6:11 PM, Nash, Anthony wrote:
 Dear Gromacs community,

 Using enforced rotation potentials I have generated a really smooth rotation
 of my ligand the catalytic binding domain of my protein. I then put together
 a simple perl script to calculate the dihedral angle of four particular atoms
 at each time step using g_angle. Taking a 6 degree interval (for a total of
 60 windows to span the 360 degree rotation) I modified the topology with a [
 dihedral_restraint ] with the respective dihedral angle for that particular
 configuation. I plan on then running each window for 10 ns to begin with
 (according to the PMF profile I can add more time). It is from this point
 where I am a little stuck.

 I have gone through a number of posts going back to 2011 with regards to
 using umbrella sampling along a dihedral angle reaction coordinate. The
 precise post was from Justin:

 You'll have to build various configurations that correspond to different
 dihedral angles (which form the sampling windows), then restrain them.

 The energy attributed to the restraints is then stored in the .edr file. From
 these energies, you should be able to construct the energy curve over the
 sampling windows. There are examples of this in the literature, so I suspect
 you should be able to find some demonstrations of how it's applied.

 I was wondering whether gromacs now support this procedure natively? I
 noticed g_wham has an option -cycl for dihedrals but I am still not sure how
 this fits in with the traditional gromacs means of generating pull force
 output as I have only been using distance based rather than angular based
 umbrella sampling.


g_wham is still only really designed to process the output from the pull code,
not any generalized restraint.  That would be a nice features, but AFAIK it is
not in the pipeline.  You can de-bias the results like I suggested above, but
it's probably more work than it's worth.  I would suggest just using Alan
Grossfield's WHAM implementation.  There, you just need to provide the raw
numbers (in this case, the dihedral values) and the force constants used.  It
will then give you your PMF.  We have done this recently; it is very easy.

-Justin

--
==

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

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

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

==
--
Gromacs Users mailing list

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

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

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

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

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

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

[gmx-users] Umbrella sampling along a dihedral angle

[Nash, Anthony]

Dear Gromacs community,

Using enforced rotation potentials I have generated a really smooth rotation of 
my ligand the catalytic binding domain of my protein. I then put together a 
simple perl script to calculate the dihedral angle of four particular atoms at 
each time step using g_angle. Taking a 6 degree interval (for a total of 60 
windows to span the 360 degree rotation) I modified the topology with a [ 
dihedral_restraint ] with the respective dihedral angle for that particular 
configuation. I plan on then running each window for 10 ns to begin with 
(according to the PMF profile I can add more time). It is from this point where 
I am a little stuck.

I have gone through a number of posts going back to 2011 with regards to using 
umbrella sampling along a dihedral angle reaction coordinate. The precise post 
was from Justin:

You'll have to build various configurations that correspond to different 
dihedral angles (which form the sampling windows), then restrain them.

The energy attributed to the restraints is then stored in the .edr file. From 
these energies, you should be able to construct the energy curve over the 
sampling windows. There are examples of this in the literature, so I suspect 
you should be able to find some demonstrations of how it's applied.

I was wondering whether gromacs now support this procedure natively? I noticed 
g_wham has an option -cycl for dihedrals but I am still not sure how this fits 
in with the traditional gromacs means of generating pull force output as I have 
only been using distance based rather than angular based umbrella sampling. 

Many thanks to any insight.
Anthony
-- 
Gromacs Users mailing list

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

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

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

[gmx-users] Renumber atom ids and charge group in topology (.itp)

[Nash, Anthony]

Dear gromacs users,

Renumbering the atom ids in a .gro file is very straight forwards, however, 
after cutting out the first 1/2 of my protein I am having great difficulty 
aligning my .itp file upon running grompp. Till this stage, it has been far 
easier for me to remove one particular domain of my protein from my 
pre-existing .gro file than it would be to adjusting the original crystal 
structure .pdb file and go through the great pains of getting pdb2gmx to work. 
But unfortunately I am now left with the fatal error:

Atoms in the .top are not numbered consecutively from 1 (rather, atomnr = 2720, 
while at-nr = 0)

which is quite understandable.

I discovered an .awk file on the gromacs website. It is approx 5 years old and 
fails to renumber my atoms in the .itp file correctly. It fails at [ angles ]. 
I've been 'googling' for a solution but I am yet to find one.

All help is appreciated.
Thanks
Anthony

Dr Anthony Nash
Department of Chemistry
University College London
-- 
Gromacs Users mailing list

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

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

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

Re: [gmx-users] Renumber atom ids and charge group in topology (.itp)

[Nash, Anthony]

Hi Justin,

If there isn't there isn't. I was just hoping to save myself the effort of 
correcting for the missing atoms on the crystal structure and the conversion of 
HIS to HIE/HID with the amber forcefield. 

Thanks
Anthony. 

Dr Anthony Nash
Department of Chemistry
University College London

From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se 
[gromacs.org_gmx-users-boun...@maillist.sys.kth.se] on behalf of Justin Lemkul 
[jalem...@vt.edu]
Sent: 16 December 2014 15:00
To: gmx-us...@gromacs.org
Subject: Re: [gmx-users] Renumber atom ids and charge group in topology (.itp)

On 12/16/14 9:54 AM, Nash, Anthony wrote:
 Dear gromacs users,

 Renumbering the atom ids in a .gro file is very straight forwards, however, 
 after cutting out the first 1/2 of my protein I am having great difficulty 
 aligning my .itp file upon running grompp. Till this stage, it has been far 
 easier for me to remove one particular domain of my protein from my 
 pre-existing .gro file than it would be to adjusting the original crystal 
 structure .pdb file and go through the great pains of getting pdb2gmx to 
 work. But unfortunately I am now left with the fatal error:

 Atoms in the .top are not numbered consecutively from 1 (rather, atomnr = 
 2720, while at-nr = 0)

 which is quite understandable.

 I discovered an .awk file on the gromacs website. It is approx 5 years old 
 and fails to renumber my atoms in the .itp file correctly. It fails at [ 
 angles ]. I've been 'googling' for a solution but I am yet to find one.


If you're hacking off part of the protein, the most straightforward solution is
to just regenerate a new topology with pdb2gmx.

-Justin

 All help is appreciated.
 Thanks
 Anthony

 Dr Anthony Nash
 Department of Chemistry
 University College London


--
==

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

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

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

==
--
Gromacs Users mailing list

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

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

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

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

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

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

Re: [gmx-users] Renumber atom ids and charge group in topology (.itp)

[Nash, Anthony]

Running through pdb2gmx as you said. 

With regards to having the original topology files, long story short, they are 
files from a half finished project of an ex-student. Files are missing and I am 
trying to resurect the work.

Thanks
Anthony  

Dr Anthony Nash
Department of Chemistry
University College London

From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se 
[gromacs.org_gmx-users-boun...@maillist.sys.kth.se] on behalf of Justin Lemkul 
[jalem...@vt.edu]
Sent: 16 December 2014 15:18
To: gmx-us...@gromacs.org
Subject: Re: [gmx-users] Renumber atom ids and charge group in  topology
(.itp)

On 12/16/14 10:16 AM, Nash, Anthony wrote:
 Hi Justin,

 If there isn't there isn't. I was just hoping to save myself the effort of 
 correcting for the missing atoms on the crystal structure and the conversion 
 of HIS to HIE/HID with the amber forcefield.


In theory, you can certainly write a converter to subtract N atoms from every
atom number that is found in the topology, but that's a far bigger pain than
just running through pdb2gmx again, in my mind.

If you had a topology before, why are there missing atom and name issues?

-Justin

 Thanks
 Anthony.

 Dr Anthony Nash
 Department of Chemistry
 University College London
 
 From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se 
 [gromacs.org_gmx-users-boun...@maillist.sys.kth.se] on behalf of Justin 
 Lemkul [jalem...@vt.edu]
 Sent: 16 December 2014 15:00
 To: gmx-us...@gromacs.org
 Subject: Re: [gmx-users] Renumber atom ids and charge group in topology (.itp)

 On 12/16/14 9:54 AM, Nash, Anthony wrote:
 Dear gromacs users,

 Renumbering the atom ids in a .gro file is very straight forwards, however, 
 after cutting out the first 1/2 of my protein I am having great difficulty 
 aligning my .itp file upon running grompp. Till this stage, it has been far 
 easier for me to remove one particular domain of my protein from my 
 pre-existing .gro file than it would be to adjusting the original crystal 
 structure .pdb file and go through the great pains of getting pdb2gmx to 
 work. But unfortunately I am now left with the fatal error:

 Atoms in the .top are not numbered consecutively from 1 (rather, atomnr = 
 2720, while at-nr = 0)

 which is quite understandable.

 I discovered an .awk file on the gromacs website. It is approx 5 years old 
 and fails to renumber my atoms in the .itp file correctly. It fails at [ 
 angles ]. I've been 'googling' for a solution but I am yet to find one.


 If you're hacking off part of the protein, the most straightforward solution 
 is
 to just regenerate a new topology with pdb2gmx.

 -Justin

 All help is appreciated.
 Thanks
 Anthony

 Dr Anthony Nash
 Department of Chemistry
 University College London


 --
 ==

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

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

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

 ==
 --
 Gromacs Users mailing list

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

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

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


--
==

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

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

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

==
--
Gromacs Users mailing list

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

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

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

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

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

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

[gmx-users] Amber to Gromacs

[Nash, Anthony]

Hi All,

I've been thrown upon a project which requires the use of the Amber FF. I have 
a crystal structure .pdb and I wish to make a topology file using the AMBER 
ff99SB forcefield. The gromacs website directs me to the ffamber ports program, 
which seems to require Gromacs versions 3.1.4, 3.2.1, . 4.0. I, however, 
have access to Gromacs 4.6.1. Must I install an old version of Gromacs? If I 
must, what kind of issues am I going to face when I take my 3.* generated 
Gromacs topology file and run it on version 4.6.1? Or is this very out of date 
and there is something new?

Thanks for any help. Apologies if these seem quite basic or I am a bit behind 
with the times, I've been doing DFT calculations for the last year and I am 
trying to get my Gromacs hat back on!

Thanks
Anthony 
-- 
Gromacs Users mailing list

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

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

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

Re: [gmx-users] Amber to Gromacs

[Nash, Anthony]

Wow! I am out of date. 

As always, thanks for the help Justin.

From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se 
[gromacs.org_gmx-users-boun...@maillist.sys.kth.se] on behalf of Justin Lemkul 
[jalem...@vt.edu]
Sent: 30 April 2014 13:29
To: gmx-us...@gromacs.org
Subject: Re: [gmx-users] Amber to Gromacs

On 4/30/14, 8:25 AM, Nash, Anthony wrote:
 Hi All,

 I've been thrown upon a project which requires the use of the Amber FF. I 
 have a crystal structure .pdb and I wish to make a topology file using the 
 AMBER ff99SB forcefield. The gromacs website directs me to the ffamber ports 
 program, which seems to require Gromacs versions 3.1.4, 3.2.1, . 4.0. I, 
 however, have access to Gromacs 4.6.1. Must I install an old version of 
 Gromacs? If I must, what kind of issues am I going to face when I take my 3.* 
 generated Gromacs topology file and run it on version 4.6.1? Or is this very 
 out of date and there is something new?


Most of the Amber force fields are supported natively now; there is no need to
download the ffamber ports so that information is indeed outdated.  Amber99-SB
and Amber99-SB-ILDN are both built into Gromacs.

-Justin

--
==

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

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

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

==
--
Gromacs Users mailing list

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

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

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

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

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

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

Re: [gmx-users] Amber to Gromacs

[Nash, Anthony]

Hi Justin,

I wonder if you could advice, I am a little further ahead than I was before. I 
am using an crystal structure .pdb from the protein data bank. There are no 
hydrogen atoms. I assume, given the pdb2gmx that this is taken care of unless 
the -ignh is selected.

I am getting a painful error upon running pdb2gmx_d -f 4AUO_mono.pdb -water 
tip3p

Opening force field file 
/Users/Anthony/gromacs451/share/gromacs/top/amber99sb.ff/aminoacids.arn
Opening force field file 
/Users/Anthony/gromacs451/share/gromacs/top/amber99sb.ff/dna.arn
Opening force field file 
/Users/Anthony/gromacs451/share/gromacs/top/amber99sb.ff/rna.arn
Checking for duplicate atoms


Processing chain 1 'A' (2986 atoms, 367 residues)
There are 568 donors and 543 acceptors
There are 786 hydrogen bonds
Will use HISE for residue 94
Will use HISE for residue 113
Will use HISE for residue 149
Will use HISE for residue 164
Will use HISE for residue 177
Will use HISE for residue 194
Will use HISE for residue 199
Will use HISE for residue 203
Will use HISE for residue 209
Will use HISE for residue 339
Will use HISE for residue 357
Will use HISE for residue 398
Will use HISE for residue 405
Will use HISE for residue 421
Identified residue PHE81 as a starting terminus.
Identified residue CYS447 as a ending terminus.
8 out of 8 lines of specbond.dat converted successfully
Special Atom Distance matrix:
   HIS94  HIS113  MET141  HIS149  HIS164  HIS177  HIS194
  NE2117  NE2275   SD493  NE2559  NE2668  NE2754  NE2916
..
...
..  [DISTANCE MATRIX HERE]
..
..
Linking CYS-259 SG-1396 and CYS-447 SG-2985...
Opening force field file 
/Users/Anthony/gromacs451/share/gromacs/top/amber99sb.ff/aminoacids.arn
Opening force field file 
/Users/Anthony/gromacs451/share/gromacs/top/amber99sb.ff/dna.arn
Opening force field file 
/Users/Anthony/gromacs451/share/gromacs/top/amber99sb.ff/rna.arn
Checking for duplicate atoms

---
Program pdb2gmx_d, VERSION 4.5.1
Source code file: pgutil.c, line: 88

Fatal error:
Atom CG not found in residue seq.nr. 11 while adding atom
---


The first residue entry is:
ATOM  1  N   PHE A  81  17.085 132.505  30.006  1.00 16.87   N
ATOM  2  CA  PHE A  81  16.981 132.301  28.532  1.00 16.72   C
ATOM  3  C   PHE A  81  18.222 132.817  27.799  1.00 16.36   C
ATOM  4  O   PHE A  81  18.921 133.701  28.291  1.00 17.46   O
ATOM  5  CB  PHE A  81  15.724 133.005  27.999  1.00 15.39   C
ATOM  6  CG  PHE A  81  15.813 134.516  27.977  1.00 11.26   C
ATOM  7  CD1 PHE A  81  16.603 135.176  27.035  1.00  8.12   C
ATOM  8  CD2 PHE A  81  15.060 135.275  28.858  1.00 10.00   C
ATOM  9  CE1 PHE A  81  16.640 136.574  26.968  1.00  6.93   C
ATOM 10  CE2 PHE A  81  15.092 136.678  28.794  1.00 10.62   C
ATOM 11  CZ  PHE A  81  15.885 137.326  27.844  1.00  7.00   C



Yet the Amberff99sb aminoacid.rtp file describes PHE as:

[ PHE ]
 [ atoms ]
 NN   -0.41570 1
 HH0.27190 2
CACT  -0.00240 3
HAH1   0.09780 4
CBCT  -0.03430 5
   HB1HC   0.02950 6
   HB2HC   0.02950 7
CGCA   0.01180 8
   CD1CA  -0.12560 9
   HD1HA   0.1330010
   CE1CA  -0.1704011
   HE1HA   0.1430012
CZCA  -0.1072013
HZHA   0.1297014
   CE2CA  -0.1704015
   HE2HA   0.1430016
   CD2CA  -0.1256017
   HD2HA   0.1330018
 CC0.5973019
 OO   -0.5679020


I am guessing that there is a bunch of atom mismatches. I am not really sure 
where to go from here. 

Many thanks
Anthony 










From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se 
[gromacs.org_gmx-users-boun...@maillist.sys.kth.se] on behalf of Nash, Anthony 
[anthony.n...@warwick.ac.uk]
Sent: 30 April 2014 14:09
To: gmx-us...@gromacs.org
Subject: Re: [gmx-users] Amber to Gromacs

Wow! I am out of date.

As always, thanks for the help Justin.

From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se 
[gromacs.org_gmx-users-boun...@maillist.sys.kth.se] on behalf of Justin Lemkul 
[jalem...@vt.edu]
Sent: 30 April 2014 13:29
To: gmx-us...@gromacs.org
Subject: Re: [gmx-users] Amber to Gromacs

On 4/30/14, 8:25 AM, Nash, Anthony wrote:
 Hi All,

 I've been thrown upon a project which requires the use of the Amber FF. I 
 have a crystal structure .pdb and I wish to make a topology file using the 
 AMBER ff99SB forcefield. The gromacs website directs me to the ffamber ports 
 program