[gmx-users] Simulation of a solid surface within Gromacs

[Kamps, M.]

Dear GMX-users,


First of all, I'm sorry if my question is a stupid one. I'm pretty new to
Gromacs and I'm trying to do my best to figure it all out.

I'm trying to work with a solid surface within Gromacs. The final goal of
my simulation is to analyse the behaviour of certain proteins on a metallic
or ceramic substrate and to find the energies related to this system. I'm
having difficulties when trying to implement a surface in Gromacs.


When I create a surface in, for instance, Avogadro or VMD Inorganic
Builder, it is possible to create a .pdb file of such a structure. However,
none of the forcefields can handle the residues, and the error: "Residue
'XXX' not found in residue topology database" arises. The User Manual gives
several options here, but I'm not sure how to approach this problem. In the
manual these options are mentioned:


*If you want a topology
 for an
arbitrary molecule, you cannot use pdb2gmx
 (unless
you build the .rtp
 entry
yourself). You will have to build it by hand, or use another program (such
as x2top  or
one of the scripts contributed by users
) to build the .top
file .*
*- see if there is a different name being used for the residue
 in the residue
database 
and rename as appropriate,*

*- parameterize
 the residue
/ molecule yourself (lots of work, even for an expert),*

*- find a topology file
 for the
molecule, convert it to an .itp file
 and
include it in your .top file
,*

*- use another force field
 which has
parameters available for this,*


*- search the primary literature for publications for parameters for the
residue that are consistent with the force field that is being used*

Based on the first sentence, I cannot use pdb2gmx and therefore have to
build the topology myself. How can I do such a thing? Is there a (simple)
guide available?

Also, am I taking the right steps for creating a surface within Gromacs? Is
there a (easier) different way?

Any help is appreciated,


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.

Re: [gmx-users] Simulation of a solid surface within Gromacs (Justin Lemkul)

[Kamps, M.]

Dear Justin, GMX-users,

First of all, thank you for the reply.
I'm not really sure how to continue from this, there is no entry in the
manual for either the command x2top or .n2t files. Google does provide some
information but only limited.
According to the information on x2top the supported input files are:

 -f  [<.gro/.g96/...>]  (conf.gro)
   Structure file: gro g96 pdb brk ent esp tpr

Where should I use this n2t file, since it is not listed? I've found this
small tutorial:

http://chembytes.wikidot.com/grocnt , with this example n2t file:

; Oplsaa-based n2t for carbon-based structures such as CNTs and graphenes
; Andrea Minoia
HHJ0.00   1.008  1C 0.109
;Hydrogen
CCJ0.00  12.011  3C 0.142   H 0.109   H 0.109
;Periferic C
CCJ0.00  12.011  3C 0.142   C 0.142   H 0.108
;Periferic C
CCJ0.00  12.011  1C 0.142
;Internal/periodic C
CCJ0.00  12.011  2C 0.142   C 0.142
;Internal/periodic C
CCJ0.00  12.011  3C 0.142   C 0.142   C 0.142
;Internal/periodic C

According to this file, the n2t file consists only of the molecules, its
name, its charge, weight, number of connections, and connected atoms? Is
this correct?

I'm kind of stuck as to how to proceed. I'm trying to implement a surface
into gromacs. This surface is obviously inorganic and can either be
metallic or ceramic, bonded or non-bonded. Is this possible?

Mark




On 11/1/16 8:44 AM, Kamps, M. wrote:
> Dear GMX-users,
>
> First of all, I'm sorry if my question is a stupid one. I'm pretty new to
> Gromacs and I'm trying to do my best to figure it all out.
>
> I'm trying to work with a solid surface within Gromacs. The final goal of
> my simulation is to analyse the behaviour of certain proteins on a
metallic
> or ceramic substrate and to find the energies related to this system. I'm
> having difficulties when trying to implement a surface in Gromacs.
>
> When I create a surface in, for instance, Avogadro or VMD Inorganic
> Builder, it is possible to create a .pdb file of such a structure.
However,
> none of the forcefields can handle the residues, and the error: "Residue
> 'XXX' not found in residue topology database" arises. The User Manual
gives
> several options here, but I'm not sure how to approach this problem. In
the
> manual these options are mentioned:
>
> *If you want a topology for an arbitrary molecule, you cannot use pdb2gmx
> (unless you build the .rtp entry yourself). You will have to build it by
hand,
> or use another program (such as x2top or one of the scripts contributed
by
> users to build the .top file
> *- see if there is a different name being used for the residue in the
residue
>   database and rename as appropriate,*
> *- parameterize the residue / molecule yourself (lots of work, even
for an
>   expert),*
> *- find a topology file for the molecule, convert it to an .itp file
and include
>   it in your .top file,*
> *- use another force field which has parameters available for this,*
> *- search the primary literature for publications for parameters for
the
>   residue that are consistent with the force field that is being used*
>
> Based on the first sentence, I cannot use pdb2gmx and therefore have to
> build the topology myself. How can I do such a thing? Is there a (simple)
> guide available?
>
> Also, am I taking the right steps for creating a surface within Gromacs?
Is
> there a (easier) different way?
>

x2top is the only GROMACS tool that will be of any use here.  You will have
to
write a suitable .n2t file (see the manual) for whatever species you're
dealing
with.

-Justin
-- 
Gromacs Users mailing list

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

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

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

[gmx-users] No interaction between manually added atoms

[Kamps, M.]

Dear GMX-users,

I'm experiencing difficulty with my manually added atoms to create a
metallic surface within GROMACS. My goal is to create a surface (at this
point metallic, Al) and let it interact with a protein.
In my current situation I've only modelled the aluminium surface
(non-bonded surface). The EM drops to a minimum in only 8 steps (first red
flag?), the NVT and NPT equilibrium look ok, but during the MD run the
surface disassembles like it is made up of loose sand. There is no
interaction at all between my individual atoms. What did I do wrong? Or
what did I forget to do?

I created the surface by defining one single Aluminium atom in a pdb file,
which i converted via pdb2gmx. To enable this, I added a Al.rtp file to a
CHARMM27 forcefield, this file contains these lines:

[ bondedtypes ]
; bonds  angles  dihedrals  impropers
 1   1  1  1
[ AL ]
 [ atoms ]
; name  type  charge  chargegroup
ALAL   0.000 0

Furthermore, to the atomtypes.atp file I added this line:

AL26.981539 ;Aluminium mass

To the ffnonbonded.itp file I added this line:

AL1326.9815390.000A0.2617500.30

Where the last two numbers relate to the Lennard-Jones potential. Although
these values are not based on any scientific literature (at this point),
they should show some interaction?

When I put the atom into pdb2gmx this is my output Al.gro file:

ALUMINIUM
1
1AL  AL1   0.000   0.000   0.000
   0.0   0.0   0.0

This atom is then replicated with insert-molecules to form an FCC crystal,
using an xyz coordinate file. The result is a nicely stacked (small slab
of) surface consisting of 256 atoms. The topology file looks like this:

[ moleculetype ]
; Namenrexcl
AL   3

[ atoms ]
;   nr   type  resnr residue  atom   cgnr charge   mass
typeBchargeB  massB
; residue   1 AL  rtp AL   q  0.0
 1 AL  1 AL AL  1  026.9815   ;
qtot 0

Including the usual lines such as include water models, ion models etc..
The surface is then put in a box (with 1.5 nm distance to edges, in all
directions) and solvated with tip3p water molecules. Can someone help me?

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.

[gmx-users] Strange behaviour of added atomtype in simulation

[Kamps, M.]

Dear gmx-user.

I am trying to implement an an-organic ionic-bound surface into GROMACS.
I've added the residue (which is Al, aluminium, uncharged) into an rtp
file, and added this to an edited Charmm force-field. The atom is added to
atomtypes.atp and ffnonbonded.itp. The LJ potentials are found from the
literature.
The surface is created by insert-molecules the individual atoms in a
specific FCC orientation. I then simulated the surface in a vacuum to see
how it behaves. To enable this, the surface is energy minimized and NVT and
NPT equilibriated.

The problem arises when it is simulated at room temperature (300K), where
the surface will behave like a liquid, instead of a solid. This seems
strange since the critical temperature is not reached by far. I tried
testing the surface at a lower temperature (10K) and there it behaves as
expected. At 10K the surface will stay in shape, in its optimal FCC
configuration and it will slightly vibrate as expected.
If the temperature is set to 300K, the surface will immediately deform and
atoms will move all through each other and around each other. The surface
will keep its shape, due to atomistic interactions, however it will not
stay in its optimal FCC configuration, which is what should happen.

What am I doing wrong?

The used md.mdp file is given below:
title = Al Surface MD simulation
; Run parameters
integrator = md-vv ; leap-frog integrator
nsteps = 20 ; 1 * 20 = 200 ps (0.2 ns)
dt = 0.001 ; 1 fs
; Output control
nstxout   = 1000 ; save coordinates every 1.0 ps
nstvout = 1000 ; save velocities every 1.0 ps
nstenergy = 1000 ; save energies every 1.0 ps
nstlog = 1000 ; update log file every 1.0 ps
; Bond parameters
continuation   = yes ; Do not continue
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
; Neighborsearching
cutoff-scheme   = Verlet
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm)
rvdw = 1.0 ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = cut-off ; Particle Mesh Ewald for long-range
electrostatics
; Temperature coupling is on
tcoupl = berendsen; modified Berendsen thermostat
tc-grps = AL ; one coupling group
tau_t = 0.1  ; time constant, in ps
ref_t = 300  ; reference temperature, one for each group, in K
; Pressure coupling is on
pcoupl   = berendsen ; Pressure coupling on in NPT
pcoupltype = semiisotropic; scaling with different value in
z direction
tau_p = 2.0; time constant, in ps
ref_p = 0.0 0.0; reference pressure, for both directions,
in bar
compressibility = 0.0 4.5e-5; compressibility, for both
directions, in bar^-1 (4.5e-5)
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Velocity generation
gen_vel = no ; velocity generation is off
; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
-- 
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] Strange behaviour of added atomtype in simulation

[Kamps, M.]

Dear gmx-users,

Is there anybody that can say anything about my simulation as to why
it would fail?
To clarify I've added some images. The first image is taken at t=0,
right after the equilibration steps, while the second image is taken
after around 0.08ns of simulations. This simulation was done in water,
however the same happens in vacuum.
https://postimg.org/image/7byeyjxo5/
https://postimg.org/image/6nfjwndyb/
As you can see, the relative order in the first image has completely
disappeared, atoms have left the FCC structure and the atoms behave
like a liquid. Atoms are also leaving the plane, and move through the
upper or lower PBC to rejoin the surface at the other side. Again,
this behaviour only happens at T=300K but not at T=10K. I understand
that atoms behave different at lower temperatures, however this
behaviour should not happen at these temperatures?
I've uploaded my files, so if you are interested you can take a look.
The first link contains all files and all simulation results, while
the second link contains only the files.
https://ufile.io/5f217
https://ufile.io/d3622
Any help would be great!

Mark

> Message: 6
> Date: Mon, 21 Nov 2016 21:05:02 +0100
> From: "Kamps, M." 
> To: gromacs.org_gmx-users@maillist.sys.kth.se
> Subject: [gmx-users] Strange behaviour of added atomtype in simulation
> Message-ID:
> 
> Content-Type: text/plain; charset=UTF-8
>
> Dear gmx-user.
>
> I am trying to implement an an-organic ionic-bound surface into GROMACS.
> I've added the residue (which is Al, aluminium, uncharged) into an rtp
> file, and added this to an edited Charmm force-field. The atom is added to
> atomtypes.atp and ffnonbonded.itp. The LJ potentials are found from the
> literature.
> The surface is created by insert-molecules the individual atoms in a
> specific FCC orientation. I then simulated the surface in a vacuum to see
> how it behaves. To enable this, the surface is energy minimized and NVT and
> NPT equilibriated.
>
> The problem arises when it is simulated at room temperature (300K), where
> the surface will behave like a liquid, instead of a solid. This seems
> strange since the critical temperature is not reached by far. I tried
> testing the surface at a lower temperature (10K) and there it behaves as
> expected. At 10K the surface will stay in shape, in its optimal FCC
> configuration and it will slightly vibrate as expected.
> If the temperature is set to 300K, the surface will immediately deform and
> atoms will move all through each other and around each other. The surface
> will keep its shape, due to atomistic interactions, however it will not
> stay in its optimal FCC configuration, which is what should happen.
>
> What am I doing wrong?
>
> The used md.mdp file is given below:
> title = Al Surface MD simulation
> ; Run parameters
> integrator = md-vv ; leap-frog integrator
> nsteps = 20 ; 1 * 20 = 200 ps (0.2 ns)
> dt = 0.001 ; 1 fs
> ; Output control
> nstxout   = 1000 ; save coordinates every 1.0 ps
> nstvout = 1000 ; save velocities every 1.0 ps
> nstenergy = 1000 ; save energies every 1.0 ps
> nstlog = 1000 ; update log file every 1.0 ps
> ; Bond parameters
> continuation   = yes ; Do not continue
> 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
> ; Neighborsearching
> cutoff-scheme   = Verlet
> ns_type = grid ; search neighboring grid cells
> nstlist = 10 ; 20 fs, largely irrelevant with Verlet scheme
> rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm)
> rvdw = 1.0 ; short-range van der Waals cutoff (in nm)
> ; Electrostatics
> coulombtype = cut-off ; Particle Mesh Ewald for long-range
> electrostatics
> ; Temperature coupling is on
> tcoupl = berendsen; modified Berendsen thermostat
> tc-grps = AL ; one coupling group
> tau_t = 0.1  ; time constant, in ps
> ref_t = 300  ; reference temperature, one for each group, in K
> ; Pressure coupling is on
> pcoupl   = berendsen ; Pressure coupling on in NPT
> pcoupltype = semiisotropic; scaling with different value in
> z direction
> tau_p = 2.0; time constant, in ps
> ref_p = 0.0 0.0; reference pressure, for both directions,
> in bar
> compressibility = 0.0 4.5e-5; compressibility, for both
> directions, in bar^-1 (4.5e-5)
> ; Periodic boundary conditions
> pbc = xyz ; 3-D PBC
> ; Velocity generation
> gen_vel = no ; velocity generation is off
> ; Dispersion correction
> DispCorr

[gmx-users] Dihedral calculations in topology

[Kamps, M.]

Dear GMX users,

I have a few questions about something I don’t fully understand. I am
trying to work with polymers, where the basics are written by Justin
Lemkul via this link:
https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-users/2009-March/040125.html

This works up to the point where I want to use grompp, where it gives
errors regarding the Proper Dihedral types; No default Proper Dih.
types. When consulting the lines where these errors refers to, it is
caused by the three dihedrals defined in my topology:

[ dihedrals ]
;  aiajakal functc0c1
c2c3c4c5
2 1 5 8 1
1 5 811 1
5 81112 1

Now, according to the manual table 5.5, page 138 (v. 2016.1), these
three dihedrals are defined as funct=1, which corresponds with a
proper dihedral, while by default the OPLS ff only knows R-B
dihedrals, which are funct=3. Up to this point my story is correct
right?

Then, why does pdb2gmx select this dihedral type? In the above
topology [ai aj ak al] = [1 5 8 11] refers to a [opls_135 opls_136
opls_136 opls_135] dihedral, which I believe is a [CT CT CT CT]
dihedral. This dihedral is already defined by ffbonded.itp as follows:
  CT CT CT CT  3  2.92880  -1.46440   0.20920
-1.67360   0.0   0.0 ; hydrocarbon all-atom
Am I correct about this part?

Is it valid to manually change the funct=1 in the above topology to
funct=3 to use this dihedral type? And if it is valid, how can I
prevent pdb2gmx to assign a funct=1 to the topology?

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.

Re: [gmx-users] Dihedral calculations in topology

[Kamps, M.]

Dear Mark,

Thanks for your reply!
I am replicating Polyethylene, of which the [polymer.rtp] and the
[polymer.hdb] are exactly the same as described by Justin Lemkul in
this link: 
https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-users/2009-March/040125.html

So, my rtp entries are:
; Terminal PE residue (chain "begin")
;C1 is -CH3
[ EthB ]
 [ atoms ]
C1opls_135   -0.180 1
   H11opls_1400.060 1
   H12opls_1400.060 1
   H13opls_1400.0601
C2opls_136   -0.120 2
   H21opls_1400.060 2
   H22opls_1400.060 2
 [ bonds ]
C1   H11
C1   H12
C1   H13
C1   C2
C2   H21
C2   H22
C2   +C1
; Terminal PE residue (chain "end")
;C2 is -CH3
[ EthE ]
 [ atoms ]
C1opls_136   -0.120 1
   H11opls_1400.060 1
   H12opls_1400.060 1
C2opls_135   -0.180 2
   H21opls_1400.060 2
   H22opls_1400.060 2
   H23opls_1400.0602
 [ bonds ]
C1   -C2
C1   H11
C1   H12
C1   C2
C2   H21
C2   H22
C2   H23

And the used hdb entries are:
EthB  2
3 4 H1 C1 C2 +C1
2 6 H2 C2 C1 +C1
EthE  2
2 6 H1 C1 C2 -C2
3 4 H2 C2 C1 -C2

The used pdb file is:
ATOM  1  C1  EthB1   0.000   0.000   0.000
ATOM  2  C2  EthB1   1.273   0.847   0.000
ATOM  3  C1  EthE1   2.546   0.000   0.000
ATOM  4  C2  EthE1   3.818   0.847   0.000

Which is a simple C4H10 molecule, composed out of two residues, with a
total of three dihedrals. The output given by grompp is:
Setting the LD random seed to -962199660
Generated 337431 of the 337431 non-bonded parameter combinations
Generating 1-4 interactions: fudge = 0.5
Generated 337431 of the 337431 1-4 parameter combinations
ERROR 1 [file topol_polymer.itp, line 119]:
  No default Proper Dih. types
ERROR 2 [file topol_polymer.itp, line 120]:
  No default Proper Dih. types
ERROR 3 [file topol_polymer.itp, line 121]:
  No default Proper Dih. types
Where lines 119 till 121 correspond with the three dihedrals in my
topology (of the polymer).

The output of pdb2gmx is as follows:
Opening force field file ./oplsaa_polymer.ff/aminoacids.r2b
Reading C4H10.pdb...
Read 4 atoms
Analyzing pdb file
Splitting chemical chains based on TER records or chain id changing.
There are 1 chains and 0 blocks of water and 2 residues with 4 atoms

  chain  #res #atoms
  1 ' ' 2  4

All occupancy fields zero. This is probably not an X-Ray structure
Opening force field file ./oplsaa_polymer.ff/atomtypes.atp
Atomtype 823
Reading residue database... (oplsaa_polymer)
Opening force field file ./oplsaa_polymer.ff/Polymer.rtp
Reading .rtp file without '[ bondedtypes ]' directive,
Will proceed as if the entry was:
[ bondedtypes ]
; bonds  angles  dihedrals  impropers all_dihedrals nr_exclusions
HH14  remove_dih
 1   1  1  2   0   3
   1  1

Residue 4
Sorting it all out...
Opening force field file ./oplsaa_polymer.ff/aminoacids.rtp
Residue 55
Sorting it all out...
Opening force field file ./oplsaa_polymer.ff/surface.rtp
Using default: not generating all possible dihedrals
Using default: excluding 3 bonded neighbors
Using default: generating 1,4 H--H interactions
Using default: removing proper dihedrals found on the same bond as a
proper dihedral
Residue 62Warning: file does not end with a newline, last line:
PTPT   0.000 0
Residue 63
Sorting it all out...
Opening force field file ./oplsaa_polymer.ff/Polymer.hdb
Opening force field file ./oplsaa_polymer.ff/aminoacids.hdb
Opening force field file ./oplsaa_polymer.ff/aminoacids.n.tdb
Opening force field file ./oplsaa_polymer.ff/aminoacids.c.tdb
Processing chain 1 (4 atoms, 2 residues)
Warning: Starting residue EthB1 in chain not identified as Protein/RNA/DNA.
Warning: Starting residue EthE1 in chain not identified as Protein/RNA/DNA.
Problem with chain definition, or missing terminal residues.
This chain does not appear to contain a recognized chain molecule.
If this is incorrect, you can edit residuetypes.dat to modify the behavior.
8 out of 8 lines of specbond.dat converted successfully
Checking for duplicate atoms
Generating any missing hydrogen atoms and/or adding termini.
Now there are 2 residues with 14 atoms
Making bonds...
Number of bonds was 14, now 13
Generating angles, dihedrals and pairs...
Before cleaning: 27 pairs
Before cleaning: 27 dihedrals
Making cmap torsions...
There are3 dihedrals,0 impropers,   24 angles
27 pairs,   13 bonds and 0 virtual sites
Total mass 58.124 a.m.u.
Total charge -0.000 e
Writing topology
Writing coordinate file...
-- 
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

Re: [gmx-users] Dihedral calculations in topology

[Kamps, M.]

Mark,

Thanks again for the swift reply! Stupid of me to miss that message! I
noticed some other code that was generated as output of pdb2gmx, which
I do not fully understand. Can you maybe explain some of it?

Processing chain 1 (4 atoms, 2 residues)
Warning: Starting residue EthB1 in chain not identified as Protein/RNA/DNA.
Warning: Starting residue EthE1 in chain not identified as Protein/RNA/DNA.
Problem with chain definition, or missing terminal residues.
This chain does not appear to contain a recognized chain molecule.
If this is incorrect, you can edit residuetypes.dat to modify the behavior.

I understand the two warnings. To my directory I added a file
residuetypes.dat, where I specificy that (among others) EthB and EthE
(which are the ending and beginning residues of polyethylene) are
Polymer. What does the rest of this message mean? I don't understand
how I could fix this message? Or even what it wants to do, since i've
defined my own terminal residues. Is it not possible to define my own
'group'?

residuetypes.dat:
EthBPolymer
EthEPolymer

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.

[gmx-users] Strange pressure coupling behaviour, rapid expanding box

[Kamps, M.]

Dear GMX users,

I'm experiencing difficulties with my simulation.
I have created a ionic bonded, rough surface of gold atoms (3276
atoms) in the OPLS FF, which is held together with Lennard-Jones
interactions. In the same box I've placed a C4H10 (14 atoms) molecule,
although the problem also occurs with proteins. This box is then
solved using SPCE water, although I've also tried TIP4P water.

The box size is originally 5.4x5.4x4.5 nm, and is first energy
minimized and NVT equilibrated. Up to this point the system behaves
normal and I can't see any red flags. When I then want to NPT
equilibrate the system, something strange happens. The periodic box
seems to ´explode´, since it rapidly (but constantly) expands. The
final dimensions are 5.4x5.4x24 nm after 50ps.

Analyzing the outputs of the NPT equilibrium shows some strange
things: The initial pressure starts at roughly 14000 bar and rapidly
decreases to around 3000 bar at t=5ps, after which it linearly
decreases to 2000 bar. Since the box is rapidly expanding, the volume
of the system is linearly growing, starting at roughly 160nm3 and
expanding to 850nm3 at 50ps.

I've added two screenshots, one before and one after the simulation.
I've also uploaded a screenshot of my mdp file, since it is easier to
read from an image than in here (i think).
https://s23.postimg.org/yceosj0q3/Before.png
https://s23.postimg.org/9kjceaobv/After.png
https://s24.postimg.org/xrzza3mt1/NPT_mdp.png

My system uses a semiisotrpic pcoupletype, since I want different
behaviour in xy and z direction. The reference pressure has been set
to 1 atmosphere (almost 1 bar), while the compressibility in the xy
direction is that of the material, while in the z direction it is that
of water. I only want periodic boundary conditions in the xy
directions, and both the z planes are LJ 12-6 walls. At z=0 the wall
is behaving like the same gold atoms, with the corresponding LJ
potential, while the upper wall is exerting an LJ potential of a dummy
atom (so no LJ interaction). Might this cause the error? I varied a
lot with the wall-r-linpot option, but putting in an positive,
negative or zero value does not seem to make a difference.

Are there any ideas as to what I am doing wrong? Thanks in advance,

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.

Re: [gmx-users] Strange pressure coupling behaviour, rapid expanding box

[Kamps, M.]

Dear Justin,

Thanks for your reply! I have tried the simulation again but this time
with opls_140, which is the same Hydrogen atom as is being used in the
polymer. This does use small values for sigma and epsilon. However
this made no difference, the box behaved the same.

Since I'm using semiisotropic coupling, and use different values for
the xy and z compressibility, I tried altering these values. In the
current situation the xy compressibility is derived from the
isothermal compressibility of gold, which is 5.77e-7. For the z
compressibility, I used 4.5e-5, which is the compressibility of water.
My thoughts were that maybe the surface caused a lot of pressure on
the system, therefore I greatly increased the compressibility in the
xy direction, to the point is higher than that of water (5.77e-4), and
I've decreased the compressibility in the z direction to 4.5e-6. This
only caused the system to expand less rapidly, but still only expand
in the z direction. The surface atoms (placed in their optimal
position), perfectly maintain the xy box.

I'm really unsure on what is causing this error, are there any other ideas?

Also, what happens when you set nwall=1? If there is no wall there,
the atoms can pass through the upper z plane, with just as the other
periodic planes, however, where do these atoms come back to? Since the
bottom plane is defined as a wall? How does this work? I can't find
this answer in the manual.

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.

[gmx-users] Protein-surface interaction during flow

[Kamps, M.]

Dear GMX users,

I am trying to simulate a flow of protein over and along a surface
structure. To do this, I have a box which is very long in its x
direction, relatively small in its y direction and its z dimension is
limited by an artificially created surface. These proteins will move
in the positive x direction, starting at around x=0, and move towards
the surface structure. At this point the protein will either pass
along the structure, or adhere to it through LJ/Coulomb interactions.
There will also be interactions between the (many) proteins flowing.
Essentially, I’m trying to create a flow.

However, I’m unsure on how to do this. If I set the proteins as a
acc-grps, they will continue to accelerate even after they have
adhered to the surface, essentially pulling them off. Is there a way
to set acc-grps for a specific amount of time or distance? After which
their velocities will guide them towards the surface structure? Or can
I give these proteins a starting velocity after which they just flow?
Any ideas?

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.

Re: [gmx-users] Protein-surface interaction during flow

[Kamps, M.]

Dear Micholas,

Thanks for your reply and suggestions!

Your idea at point 2, how can I do this? The forces and velocities are
stored in the .edr files, while the final velocities and positions are
also stored in the output .gro file right? This means I should run a
simulation just long enough that the proteins reach a desired point,
after which I can use the final .gro  or the .edr files?
I then need to continue (extend) my simulation, while creating a new
.tpr file via the new .mdp file, effectively using this command:

gmx grompp -f new.mdp -c old.tpr -o new.tpr -t old.cpt
mdrun -s new.tpr
(http://www.gromacs.org/Documentation/How-tos/Extending_Simulations)

Will this continue my simulation with the 'previous' velocities? Since
none of these files will input any .edr or .gro files?

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.

[gmx-users] Velocity as a function of distance Z

[Kamps, M.]

Dear gmx-users,

I have several simulations regarding atoms flowing between two gold
surfaces. The flow of the liquid molecules show a distinct Poiseuille
flow, which was expected.
I am now interested in extracting the velocities of these liquid
molecules (total of around 3 atoms, around 700 molecules).

I am interested in the (average) velocity, as a function of the height
Z. Is this possible in any way?

If I take for example the command gmx traj ... -av, I will get the
average velocity of all atoms, however I cannot tell which atom is
placed where. Furthermore, is it possible to get the average
velocities per molecule instead of atoms, for instance via gmx traj?

Thanks in advance.

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.

Re: [gmx-users] Velocity as a function of distance Z

[Kamps, M.]

Dear Mark,

Thanks for your reply. Correct me if i'm wrong, but your suggestion
would be to determine the position of each atom/molecule during a
frame. This process should be repeated for each frame, which is
easiest to do with a manually written (non-GROMACS) script.
With these positions per time, gmx traj should be able to determine
the velocity of each atom/molecule? If so, I will try to figure this
out.

Furthermore, is there a way to automatically select multiple atoms as
a molecule? My simulation uses small simple polymers, is there a way
to let GROMACS calculate only the positions or velocities of the
molecules instead of the atoms?

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.

Re: [gmx-users] Velocity as a function of distance Z (Mark Abraham)

[Kamps, M.]

Dear Mark,

Thanks again for your reply. I'm sorry for asking these probably
stupid questions, but I'm not able to figure it out.

The objects I am interested in are small polymers, consisting of three
residue groups; Eth EthE and EthB (corresponding to both the end
groups of a polymer, and the middle section). I total these groups
consist of around 25000 atoms.

With gmx select I will enter the following command, at an arbitrarily
chosen time:
gmx select -f input.trr -s input.tpr -on output.ndx -b 7.5 -e 7.5
In the following selection screen my preferred molecules are listed as
three different groups (Eth, EthE and EthB), therefore I manually
select these groups with the command: resname Eth EthE EthB, press
enter and end with Ctrl-D. It will then proceed to process the frames,
where GROMACS tells me it has analyzed 1 frame, at the 30th timestep,
which is 7.5.

I will then switch to gmx traj where I enter the following:
gmx traj -f input.trr -s input.tpr -n output.ndx -ov output.xvg
It will then process all frames, and not only the frame I selected
during gmx select. Do I need to specify -b and -e again during gmx
traj?

Then, the output could be read with xmgrace, however looking at the
files it will try to plot 25000 different lines, which it probably
cant. Opening the .xvg file plots the velocity/time plot of only 1
atom, while opening the  .xvg file for the time selected output (via
gmx traj -b -e) shows an empty plot for only one atom.

I'm getting lost in all the options, any help would be appreciated.

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.

Re: [gmx-users] gromacs.org_gmx-users Digest, Vol 155, Issue 21

[Kamps, M.]

Dear Mark,

Again, thanks for the reply.

You said I should identify groups of molecules that are to my
interest. I'm unsure how to do this.
My .gro files show no distinction between different molecules, all of
them are built the same (rows deleted for clarification).

1EthBC1 5963  18.910   0.037   2.874  0.4243  0.6722 -0.0361
1EthB   H11 5964  18.969   0.099   2.810  1.4488  0.2545  0.0251
  ...
5EthH21 5992  18.287   0.495   3.423  0.1969  0.2998  0.4129
5EthH22 5993  18.291   0.508   3.251  0.7485 -2.7065  0.0733
  ...
   10EthE   H22 6023  17.188   0.735   3.427 -0.4008 -2.3617  0.5620
   10EthE   H23 6024  17.153   0.699   3.263  2.2014  0.6628 -1.8300

What should I enter to select individual molecules? There are no
columns or names that are unique to a single molecule?
In the 'Selection syntax and usage' document there are tons of
options, however I find it difficult to know which one to use, and how
to use them.
The distinctions that I can make are in their residue
(EthB-(Eth)n-EthE), which will select all atoms in these groups,
instead of whole molecules. I could also make a selection based on
their position, however that would mean not entire molecules are
taken. When i say for example resname EthB EthE Eth and x>5 and x<7,
all atoms inside this margin will be taken, however this would mean it
will select partial/split/cut molecules, since parts of the molecule
will be outside of this frame. When looking at the earlier mentioned
document, I see for example the 'chain' option. But I can't figure out
how to use this.

My residue groups are defined via residuetypes.dat to be:





>> Dear Mark,
>>
>> Thanks again for your reply. I'm sorry for asking these probably
>> stupid questions, but I'm not able to figure it out.
>>
>> The objects I am interested in are small polymers, consisting of three
>> residue groups; Eth EthE and EthB (corresponding to both the end
>> groups of a polymer, and the middle section). I total these groups
>> consist of around 25000 atoms.
>>
>
> Those names won't help you. You want to identify (groups of) single
> molecules. Look at your file and see what distinguishes them.
>
>> With gmx select I will enter the following command, at an arbitrarily
>> chosen time:
>> gmx select -f input.trr -s input.tpr -on output.ndx -b 7.5 -e 7.5
>> In the following selection screen my preferred molecules are listed as
>> three different groups (Eth, EthE and EthB), therefore I manually
>> select these groups with the command: resname Eth EthE EthB, press
>> enter and end with Ctrl-D. It will then proceed to process the frames,
>> where GROMACS tells me it has analyzed 1 frame, at the 30th timestep,
>> which is 7.5.
>>
>
> I can't see everything, but probably all you did was select every atom,
> which won't help you. You wanted all the molecules with some criterion, so
> you will need to make a better selection.
>
>> I will then switch to gmx traj where I enter the following:
>> gmx traj -f input.trr -s input.tpr -n output.ndx -ov output.xvg
>> It will then process all frames, and not only the frame I selected
>> during gmx select. Do I need to specify -b and -e again during gmx
>> traj?
>>
>
> Yes, once you've made a selection that matches geometric criteria from a
> frame, it is literally only applicable to that frame.
>
>> Then, the output could be read with xmgrace, however looking at the
>> files it will try to plot 25000 different lines, which it probably
>> cant. Opening the .xvg file plots the velocity/time plot of only 1
>> atom, while opening the  .xvg file for the time selected output (via
>> gmx traj -b -e) shows an empty plot for only one atom.
>>
>
> You haven't selected eg one molecule from one frame yet, so get that right
> first.
>
> Mark
>
>> I'm getting lost in all the options, any help would be appreciated.
>>
>> 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.

Re: [gmx-users] Velocity as a function of distance Z

[Kamps, M.]

(previous mail was sent prematurely, my apologies)

Dear Mark,

Again, thanks for the reply.

You said I should identify groups of molecules that are to my
interest. I'm unsure how to do this.
My .gro files show no distinction between different molecules, all of
them are built the same (rows deleted for clarification).

1EthBC1 5963  18.910   0.037   2.874  0.4243  0.6722 -0.0361
1EthB   H11 5964  18.969   0.099   2.810  1.4488  0.2545  0.0251
  ...
5EthH21 5992  18.287   0.495   3.423  0.1969  0.2998  0.4129
5EthH22 5993  18.291   0.508   3.251  0.7485 -2.7065  0.0733
  ...
   10EthE   H22 6023  17.188   0.735   3.427 -0.4008 -2.3617  0.5620
   10EthE   H23 6024  17.153   0.699   3.263  2.2014  0.6628 -1.8300

What should I enter to select individual molecules? There are no
columns or names that are unique to a single molecule?
In the 'Selection syntax and usage' document there are tons of
options, however I find it difficult to know which one to use, and how
to use them.
The distinctions that I can make are in their residue
(EthB-(Eth)n-EthE), which will select all atoms in these groups,
instead of whole molecules. I could also make a selection based on
their position, however that would mean not entire molecules are
taken. When i say for example resname EthB EthE Eth and x>5 and x<7,
all atoms inside this margin will be taken, however this would mean it
will select partial/split/cut molecules, since parts of the molecule
will be outside of this frame. When looking at the earlier mentioned
document, I see for example the 'chain' option. But I can't figure out
how to use this.

My residue groups are defined via residuetypes.dat to be:
EthPolymer
EthBPolymer
EthEPolymer
While in my topology they are called simply:
C20H42
C16H34
Is this helpful to make a molecule selection?


Mark




>>> Dear Mark,
>>>
>>> Thanks again for your reply. I'm sorry for asking these probably
>>> stupid questions, but I'm not able to figure it out.
>>>
>>> The objects I am interested in are small polymers, consisting of three
>>> residue groups; Eth EthE and EthB (corresponding to both the end
>>> groups of a polymer, and the middle section). I total these groups
>>> consist of around 25000 atoms.
>>>
>>
>> Those names won't help you. You want to identify (groups of) single
>> molecules. Look at your file and see what distinguishes them.
>>
>>> With gmx select I will enter the following command, at an arbitrarily
>>> chosen time:
>>> gmx select -f input.trr -s input.tpr -on output.ndx -b 7.5 -e 7.5
>>> In the following selection screen my preferred molecules are listed as
>>> three different groups (Eth, EthE and EthB), therefore I manually
>>> select these groups with the command: resname Eth EthE EthB, press
>>> enter and end with Ctrl-D. It will then proceed to process the frames,
>>> where GROMACS tells me it has analyzed 1 frame, at the 30th timestep,
>>> which is 7.5.
>>>
>>
>> I can't see everything, but probably all you did was select every atom,
>> which won't help you. You wanted all the molecules with some criterion, so
>> you will need to make a better selection.
>>
>>> I will then switch to gmx traj where I enter the following:
>>> gmx traj -f input.trr -s input.tpr -n output.ndx -ov output.xvg
>>> It will then process all frames, and not only the frame I selected
>>> during gmx select. Do I need to specify -b and -e again during gmx
>>> traj?
>>>
>>
>> Yes, once you've made a selection that matches geometric criteria from a
>> frame, it is literally only applicable to that frame.
>>
>>> Then, the output could be read with xmgrace, however looking at the
>>> files it will try to plot 25000 different lines, which it probably
>>> cant. Opening the .xvg file plots the velocity/time plot of only 1
>>> atom, while opening the  .xvg file for the time selected output (via
>>> gmx traj -b -e) shows an empty plot for only one atom.
>>>
>>
>> You haven't selected eg one molecule from one frame yet, so get that right
>> first.
>>
>> Mark
>>
>>> I'm getting lost in all the options, any help would be appreciated.
>>>
>>> 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.

[gmx-users] Meaning of columns in .xvg output file

[Kamps, M.]

Dear gmx-users,

I've got a small question of which I can't find the answer. I have a
simulation where I want to calculate the average velocity of certain
particles. I've created an output .xvg file via the command:
gmx traj -f input.trr -s input.tpr -n input.ndx -av output.xvg
My output, is as requested an .xvg file which can be plotted in xmgrace.

However, when opening this .xvg file with a text-editor, it has four
columns of output, of which only the first two columns will be used by
xmgrace. The first column is a simple counter, however the other three
columns all contain values. Supposedly, the second column (column 2 of
4) contains the average velocities. My question: What are the values
in column 3 and 4?

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.

Re: [gmx-users] Meaning of columns in .xvg output file

[Kamps, M.]

Dear Justin,

Thank you for you reply. The header of my .xvg file is:

# This file was created Mon Mar  6 14:32:28 2017
# Created by:
#   :-) GROMACS - gmx traj, 2016.1 (-:
#
# Executable:   /usr/local/gromacs/bin/gmx
# Data prefix:  /usr/local/gromacs
# Working dir:  /local/Data/Nozzle/19./4.MD Acceleration
# Command line:
#   gmx traj -f md_acc.trr -s md_acc.tpr -b 7.5 -e 8.5 -n 1.ndx -av plot1.xvg
# gmx traj is part of G R O M A C S:
#
# GROtesk MACabre and Sinister
#
@title "average velocity"
@xaxis  label "Atom"
@yaxis  label ""
@TYPE xy
1   0.015   0.354  -0.083
2  -0.242  -0.053   0.026

Where the counter obviously continues to larger values. As you can see
the x-axis contains the atom number (same as the first column; the
counter), while the second column is plotted as the y-values. The
y-axis remains empty.
The .xyz file is created via gmx select to select multiple residue
groups, after which I am interested in the average velocity of the
selected groups. Is it possible that the three colums represent the
average velocity in each of the three directions (XYZ?)

Mark





> On 3/6/17 5:36 AM, Kamps, M. wrote:
>> Dear gmx-users,
>>
>> I've got a small question of which I can't find the answer. I have a
>> simulation where I want to calculate the average velocity of certain
>> particles. I've created an output .xvg file via the command:
>> gmx traj -f input.trr -s input.tpr -n input.ndx -av output.xvg
>> My output, is as requested an .xvg file which can be plotted in xmgrace.
>>
>> However, when opening this .xvg file with a text-editor, it has four
>> columns of output, of which only the first two columns will be used by
>> xmgrace. The first column is a simple counter, however the other three
>
> You can plot everything with xmgrace -nxy output.xvg
>
>> columns all contain values. Supposedly, the second column (column 2 of
>> 4) contains the average velocities. My question: What are the values
>> in column 3 and 4?
>>
>
> What do the labels in the header of the file say, or alternatively, what does
> the legend say when you plot with the above command?
>
> -Justin
-- 
Gromacs Users mailing list

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

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

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

[gmx-users] acc-grps and acceleration

[Kamps, M.]

Dear GMX-users,

I have a question relating the acceleration of certain atoms in my system.
I am trying to simulate a flow of a liquid between two gold plates.
I've created a small, narrow and long box filled with this liquid.

In order to create the flow, a acceleration is applied to the liquid.
However I'm having difficulty understanding the behaviour of this
acceleration. I understand that, due to computational reasons, the
applied acceleration is very high (0.1 nm/ps^2 = 1x10^14 m/s^2).
However, due to the small time-scale, this acceleration still leads to
usable results (literature indicates this).

I'm having difficulty with the development of the flow. I would
suspect to see the flow velocity develop over time. I would suspect
that at t=0 the average velocity is zero, while at t=a the flow would
have a velocity which is non-zero. At t=b the applied acceleration
(=constant force) would be overcome by the fluid/wall friction, and a
steady-state flow would emerge. However, the behaviour that I'm seeing
does not show this 'development' in flow velocity. At t=0 the velocity
is already at a higher velocity and this velocity does not increase
over time.

I created my system and thoroughly applied EM and equilibrated it. In
my final equilibration step I generated a velocity through gen_vel,
while no acceleration was applied. For my final MD run no velocity was
generated and the final simulation was continued from this
equilibration step, with an acceleration.

I've tried to apply different velocities (1 0.1 0.01 0.001, all in
correct direction), and I've varied with the length of the simulation
(500ps, 50ps, 5ps and 0.5ps) and with the dt value (0.001 0.0005 and
even 0.5).

Any ideas as to what is happening? Or is this behaviour normal and my
interpretation wrong?
-- 
Gromacs Users mailing list

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

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

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

[gmx-users] Filling box with long linear molecules

[Kamps, M.]

Dear GMX users,

I am working with polymers (~20 carbon atoms) in an enclosed box. In
order to fill this box with polymers, I currently use an iteration of
multiple steps.

Currently, the polymers are created via an external code, and the
resulting polymers are linear strings of polymers, with only a
'zig-zag' motion at the carbon atoms. Currently I am filling this box
via the 'gmx insert-molecules -nmol X' command, where it will try to
fill the box.

Due to the orientation of the polymers, the box is considered full and
GROMACS cannot add any more polymers. However, after performing a
small relaxation the polymers will 'bend and twist' all around each
other, resulting in lots of empty space. This process is repeated
several times until gmx insert-molecules can no longer fit any more
polymers in the box.

My question, is there a more efficient way to do this? Can GROMACS
insert-molecules while 'bending' the polymers?

Kind regards,

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.

[gmx-users] Unreasonably high pressure values

[Kamps, M.]

Dear GMX users,

I am struggling with understanding my systems behaviour.
I have a semi-isotropic system where my bottom and upper plane are
composed of gold FCC plates. In between is a fluid, which during
equilibration is stationary (no acceleration). I equilibrate the
system with a Berendsen barostat, as the literature suggests. The
compressibility is that of the fluid used (Hexadecane C16H34), while
the reference pressure is 1 bar. See the below section of my systems
MDP:

; Pressure coupling is on
pcoupl= berendsen; Pressure coupling on
pcoupltype= semiisotropic; uniform scaling
of box vectors
tau_p= 2.0; time constant, in ps
ref_p= 1.0 1.0; reference pressure, in bar
compressibility = 0  8.6e-5; isothermal
compressibility, bar^-1

When I equilibrate the system, everything appears to be normal. The
volume shrinks a little to compress the fluid, which leads to a
correct density (correct as in corresponding with known values of the
real-life density). The fluid behaves as expected and everything looks
good.

I then want to use the system to accelerate the fluid, however that
blows up the system. Blowing up as in; fluctuations in the box size,
shrink, expand, rapid shrink, rapid expand, extreme shrink, blowing
up.

Upon inspection, during equilibrating the pressure of my system nicely
converges, however, it converges to -6000 bar, which seems rather
strange.. How does this happen? What am I doing wrong?

I use a 0.5 fs timestep, equilibriate over 400.000 steps (0.2 ns),
300K v-rescale temperature coupling, PME electrostatic calculations.

Kind regards,
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.

Re: [gmx-users] Unreasonably high pressure values

[Kamps, M.]

Dear Mark,

Thanks for your answer. I did some experiments and I think I've
located the problem.
Within my simulation there are two main components: Two surfaces and a
polymer fluid. These two surfaces are created by placing gold atoms in
a FCC lattice (lattice parameters derived from the LJ potential). The
only interactions between these surfaces are the LJ potentials.

When I run a simple NVT equilibrium on a basic system (without
pressure, without accelerations, etc), the pressure is already really
high (again, around 6000 bar, it stabilizes around this value!). When
I run the same NVT equilibrium on only the fluid (so no surfaces), the
pressure remains stable at around 0-1 bar, as should be expected.

Performing an NPT equilibrium after this first one grants the same
results. The system with the surfaces shows unreasonably high pressure
values, while the fluid-only system behaves as expected. (The
pressures are calculated via gmx energy)

My guess, the problem is related to the gold surfaces. However, when
running only a NVT and NPT equilibration on the surfaces, nothing
happens. There are no high stresses within the system, everything runs
fine. The only difference is the pressure, which is extremely high.

Is there a way to exclude the surface from the calculations? I am only
interested in the LJ behaviour of the surface towards the fluid. I
cannot use the GROMACS walls options since I am interested in non-flat
surface.

Any ideas?

___

Message: 1
Date: Mon, 15 May 2017 16:35:56 +
From: Mark Abraham 
To: gmx-us...@gromacs.org, gromacs.org_gmx-users@maillist.sys.kth.se
Subject: Re: [gmx-users] Unreasonably high pressure values
Message-ID:
 
Content-Type: text/plain; charset="UTF-8"

Hi,

Convergence times depend on system type and size, of course, but an
equilibration on an inhomogeneous system should probably run for more than
a few nanoseconds. That might help expose whether the issue is the system,
or the protocol, or the applied force.

Mark

On Mon, 15 May 2017 18:25 Kamps, M.  wrote:

> Dear GMX users,
>
> I am struggling with understanding my systems behaviour.
> I have a semi-isotropic system where my bottom and upper plane are
> composed of gold FCC plates. In between is a fluid, which during
> equilibration is stationary (no acceleration). I equilibrate the
> system with a Berendsen barostat, as the literature suggests. The
> compressibility is that of the fluid used (Hexadecane C16H34), while
> the reference pressure is 1 bar. See the below section of my systems
> MDP:
>
> ; Pressure coupling is on
> pcoupl= berendsen; Pressure coupling on
> pcoupltype= semiisotropic; uniform scaling
> of box vectors
> tau_p= 2.0; time constant, in ps
> ref_p= 1.0 1.0; reference pressure, in
> bar
> compressibility = 0  8.6e-5; isothermal
> compressibility, bar^-1
>
> When I equilibrate the system, everything appears to be normal. The
> volume shrinks a little to compress the fluid, which leads to a
> correct density (correct as in corresponding with known values of the
> real-life density). The fluid behaves as expected and everything looks
> good.
>
> I then want to use the system to accelerate the fluid, however that
> blows up the system. Blowing up as in; fluctuations in the box size,
> shrink, expand, rapid shrink, rapid expand, extreme shrink, blowing
> up.
>
> Upon inspection, during equilibrating the pressure of my system nicely
> converges, however, it converges to -6000 bar, which seems rather
> strange.. How does this happen? What am I doing wrong?
>
> I use a 0.5 fs timestep, equilibriate over 400.000 steps (0.2 ns),
> 300K v-rescale temperature coupling, PME electrostatic calculations.
>
> Kind regards,
> 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.

Re: [gmx-users] Unreasonably high pressure values

[Kamps, M.]

Dear Mark,

I'm sorry for being unclear. I have rectangular box (20x5x7 nm, XYZ)
where a flow should be replicated. at Z=0 and Z=7 there are two
surfaces composed of gold atoms. In between there is a fluid. The
system is open or periodic in the X and Y directions, but limited in
the Z direction by the surfaces (PBC is set to XYZ, while periodic
molecules is set to yes). Both surfaces 'connect' through the PBC. The
goal is to replicate a fluid flow between the two surfaces, which will
by done by applying a constant acceleration to the fluid in the X
direction. The flow will be therefore along the X axis. It is somewhat
similar to articles as these:
http://advances.sciencemag.org/content/2/3/e1501585/tab-pdf

I don't fully understand what you mean in your last sentences. Why
shouldn't my system produce sane values of pressure?


___


Message: 1
Date: Tue, 16 May 2017 09:11:46 +
From: Mark Abraham 
To: gmx-us...@gromacs.org, gromacs.org_gmx-users@maillist.sys.kth.se
Subject: Re: [gmx-users] Unreasonably high pressure values
Message-ID:

Content-Type: text/plain; charset="UTF-8"

Hi,

I don't understand the system very well from your description. You have two
surfaces, a fluid between them, presumably periodic in the dimensions of
the surface. But what do you have in the direction normal to the surface?
Whatever it is, why should your parameterization of your gold atoms LJ lead
to sane pressure? What do you get from a normal 3D periodic system with a
surface solvated on both sides?

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.

[gmx-users] Problem with accelerations

[Kamps, M.]

Dear GMX users,

I have some strange behaviour which I cannot explain.

I want to accelerate atoms through my box at a certain velocity. Since
I can only adjust the acceleration, I have to trial-and-error my way
to the right accelerations.

To do this, I create a smaller 'testing' simulation, which is
continued from an extensive equilibrium. I apply an acceleration of
0.05 nm/ps2, and after 200ps my velocity (gathered from gmx traj with
an gmx select input) says that the velocity is stable at around 0.04
nm/ps. throughout the simulation the velocity slightly increases due
to the atoms rearranging etc.

Now I want to simulate the same behaviour, but for a longer amount of
time. I therefore take the exact same MDP file, and change nothing
except the time-related parameters. I change the number of total steps
and the timestep, but leave the acceleration intact. After analysing
the data I get a MUCH higher velocity! I can understand this due to
rearranging of the atoms on the longer term, but after the same 200ps
the velocity is also way higher.

So during my smaller run the velocity was 0.04 nm/ps after 200 ps, in
this longer run, the velocity is 0.2 nm/ps after the same 200ps. How
is this possible? The exact same acceleration is applied.

Am I missing something?
-- 
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] Problem with accelerations

[Kamps, M.]

Mark,

I'm not really sure what you mean by technique. I assume NEMD stands
for Non-equilibrium MD? That is the case..

About my simulation: I am trying to simulate a fluid flow between
surfaces. Two atomistic surfaces of gold atoms have been created at
Z=0 and Z=z, where z is the top of the box. These surfaces are fixed.
In between is a fluid that should be accelerated in order to create a
flow. This is done by applying a constant acceleration in the X
direction via the commands; 'acc-grps = FLUID' and 'accelerate = 0.05
0 0'. Is this equal to the bug you referenced to?

Suppose the command is broken, does this mean that the simulation is
worthless since there are computational errors, or is the applied
velocity buggy but the results still useful? Suppose I can find the
right accelerations, can I still use the results?


Mark Abraham wrote:
> Hi,
> Which technique are you using for this?
> https://redmine.gromacs.org/issues/1354 speculates that this code is just
> broken, has been so for years, and should be removed given that nobody
> wishes to do the work to fix it ( :-( )
> Mark
> On Tue, May 23, 2017 at 11:13 AM Kamps, M.  wrote:
> > Dear GMX users,
> >
> > I have some strange behaviour which I cannot explain.
> >
> > I want to accelerate atoms through my box at a certain velocity. Since
> > I can only adjust the acceleration, I have to trial-and-error my way
> > to the right accelerations.
> >
> > To do this, I create a smaller 'testing' simulation, which is
> > continued from an extensive equilibrium. I apply an acceleration of
> > 0.05 nm/ps2, and after 200ps my velocity (gathered from gmx traj with
> > an gmx select input) says that the velocity is stable at around 0.04
> > nm/ps. throughout the simulation the velocity slightly increases due
> > to the atoms rearranging etc.
> >
> > Now I want to simulate the same behaviour, but for a longer amount of
> > time. I therefore take the exact same MDP file, and change nothing
> > except the time-related parameters. I change the number of total steps
> > and the timestep, but leave the acceleration intact. After analysing
> > the data I get a MUCH higher velocity! I can understand this due to
> > rearranging of the atoms on the longer term, but after the same 200ps
> > the velocity is also way higher.
> >
> > So during my smaller run the velocity was 0.04 nm/ps after 200 ps, in
> > this longer run, the velocity is 0.2 nm/ps after the same 200ps. How
> > is this possible? The exact same acceleration is applied.
> >
> > Am I missing something?
> > --
> > 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] Problem with accelerations

[Kamps, M.]

Mark,

I'm not really sure what you mean by technique. I assume NEMD stands
for Non-equilibrium MD? That is the case in my simulation..

About my simulation: I am trying to simulate a fluid flow between
surfaces. Two atomistic surfaces of gold atoms have been created at
Z=0 and Z=z, where z is the top of the box. These surfaces are fixed.
In between is a fluid that should be accelerated in order to create a
flow. This is done by applying a constant acceleration in the X
direction via the commands; 'acc-grps = FLUID' and 'accelerate = 0.05
0 0'. Is this equal to the bug you referenced to?

Suppose the command is broken, does this mean that the simulation is
worthless since there are computational errors, or is the applied
velocity buggy but the results still useful? Suppose I can find the
right accelerations, can I still use the results?

Another question, how can I properly 'reply' to an answer working from
the digest mailing list? If I reply normally (via my mail client) and
delete the other messages, my message will not be displayed as a reply
in online mailing-list viewers such as Narkive etc. Am I replying
correctly?


Mark Abraham wrote:
> Hi,
> Which technique are you using for this?
> https://redmine.gromacs.org/issues/1354 speculates that this code is just
> broken, has been so for years, and should be removed given that nobody
> wishes to do the work to fix it ( :-( )
> Mark
> On Tue, May 23, 2017 at 11:13 AM Kamps, M.  wrote:
> > Dear GMX users,
> >
> > I have some strange behaviour which I cannot explain.
> >
> > I want to accelerate atoms through my box at a certain velocity. Since
> > I can only adjust the acceleration, I have to trial-and-error my way
> > to the right accelerations.
> >
> > To do this, I create a smaller 'testing' simulation, which is
> > continued from an extensive equilibrium. I apply an acceleration of
> > 0.05 nm/ps2, and after 200ps my velocity (gathered from gmx traj with
> > an gmx select input) says that the velocity is stable at around 0.04
> > nm/ps. throughout the simulation the velocity slightly increases due
> > to the atoms rearranging etc.
> >
> > Now I want to simulate the same behaviour, but for a longer amount of
> > time. I therefore take the exact same MDP file, and change nothing
> > except the time-related parameters. I change the number of total steps
> > and the timestep, but leave the acceleration intact. After analysing
> > the data I get a MUCH higher velocity! I can understand this due to
> > rearranging of the atoms on the longer term, but after the same 200ps
> > the velocity is also way higher.
> >
> > So during my smaller run the velocity was 0.04 nm/ps after 200 ps, in
> > this longer run, the velocity is 0.2 nm/ps after the same 200ps. How
> > is this possible? The exact same acceleration is applied.
> >
> > Am I missing something?
> > --
> > 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] Strange acceleration behaviour after continuation

[Kamps, M.]

Dear GMX users,

I am trying to replicate a molecular flow of small molecules between
two solid slabs of FCC bound atoms. The slabs are bound via LJ
interactions.The flow is created by using periodicity, where the fluid
is accelerated continuously between the two slabs. The periodic box is
created as a long small 'channel', where periodicity can occur within
the X and Y directions. The Z direction is restricted due to the fixed
'walls'.

I am now trying to control the fluid flow, by altering the
accelerations. To do this, I first accelerate the fluid with a
constant acceleration 'a1'. This will result in a movement of the
fluid over time. At t=0 the velocity is zero, while at t=1 the
velocity equals 0.06 nm/ps2. The velocity development is linear,
meaning it has a linear increase from 0 to this value of 0.06 nm/ps2.

After this 'first stage', I want a less rigorous increase within the
velocity, therefore I perform a second simulation, continued from the
previous simulation. In this case the acceleration is set to a value
that is 20% that of 'a1' (the previous acceleration). Within the MDP
file it is set to genvel=no and continuation=yes. The .cpt and .gro
files are provided as the input for this second simulation.

My expectation is that the velocities at the beginning of this
simulation are equal to that of previous simulation, hence the
continuation. However, this is not the case. The velocities at the
beginning of this simulation start at roughly 0.02 nm/ps2 and increase
(slowly) to 0.04 nm/ps.

My question: why is there a mismatch in the velocities? Am I doing
something wrong? Is my logic flawed?

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.