You might try building the initial structure and
defining the force field in a different molecular
simulation package (i.e. CHARMM, AMBER). Choose a
package that has the force field you have chosen
(really you should start by picking the force field
you want to use (i.e. CHARMM, AMBER).
2. Energy Minimization to get a realistic structure
LAMMPS doesn't do this, but there are ways around
this.
3. (NPT) dynamics at high pressure for 100 ps to
arrive at a realistic
density.
2) Alternatively, you could gradually increase the
timestep size as in this example:
units real
neigh_modify delay 5 every 1
atom_style full
bond_style harmonic
angle_style harmonic
dihedral_style charmm
pair_style lj/charmm/coul/long 8 10
kspace_style pppm 1e-4
read_data my.data
special_bonds 0.0 0.0 1.0
thermo_style multi
thermo 1
velocity all create 300.0 4928459
fix 1 all temp/rescale 298.0 298.0 1 100.
0.5
fix 2 all nve
timestep 0.01
run 100
timestep 0.1
run 100
fix 3 all shake 0.000001 500 0 m 1.0 a 4
timestep 2.0
thermo_style one
thermo 500
run 500000
I'm trying to implement Siewert Marrink's coarse grain model (S.J.
Marrink, A.H. de Vries, A.E. Mark. Coarse grained model for
semi-quantitative lipid simulations. JPC-B, 108:750-760, 2004.) in LAMMPS.
I've gotten all the potential parameters ported over, but I'm having
problems with using the same timestep that he does (50 fs). I minimize
using a timestep of 1fs, 0 velocity and NVE, and then heat up
(temp/rescale) with a timestep of 10fs. Everything is fine at this point,
but when I switch to NPT with an integration timestep of 50fs I get nan's
after the 0th timestep. If I run NPT at a timestep of 10fs it works fine,
but as soon as I switch it give's me nan's. Is there any debugging output
I can print out to help me narrow down the problem?
Here are some suggestions that I hope you'll find
useful:
1) You might need to equilibrate longer with the
shorter timestep.
2) There may be something wrong with your
implementation of the force field.
3) Is it possible that they are mapping to 50 fs
timesteps, but really running 10 fs timesteps?
4) Have you tried it with NVE and 50 fs timestep?
Probably best if you can send your input files so that
I might more easily spot things that may be amiss.
Unfortunately, I don't have suggestions on how to get
LAMMPS to self-diagnose what the problem is. I don't
know of an easy way to get it to dump debugging output
to help you narrow things down. Staring at the input
files, where the problem will most likely be found,
seems to be the most effective approach to solving
problems like these.
Is there a limitation on the size of the domain that can be simulated using
LAMMPS? I'm running the program in:
Intel EM64T processor (64-bit) - dual proc; 4GB memory Linux Redhat AS 3.0
Intel compilers 8.1 HP MPI
When I try to represent a large number of atoms (box size 0 3000 0 1000 0
3000) the code seems to run without generating any outputs, e.g. I'm
interested in obtaining mean square displacements of a certain group of
atoms in the domain, but for large box sizes the program doesn't produce the
desired results.
assuming that you have defined the region without the 'units box '-parameter, it means that you are trying to generate 3000x1000x3000x(# of base atoms) which is at least 9x10^9 atoms. Your computer will be busy a long time just to generate all these atoms, and it will crash eventually as you do no have the memory (more than 700 GB ) for so many atoms.
you can also dump out forces, velocities, coords in a dump
file every timestep ... so you should quickly see what
forces are causing the system to blow up
Is it possible to specify that some bond types which are not excluded from
pairwise interactions and some that are? I am trying to implement a
harmonic distance based restraint (to model a hydrogen bond interaction)
so I would like them to not be excluded from pairwise interactions, while
excluding the actual "physically" bonded atoms. I can't use fix spring
because 1) i need a whole bunch of these and there is a limit of 32 groups
2) i actually want to model a restrain to a distance region, not a fixed
value (these are often used when refining NMR models using NOE data).
Would it be possible in the current code or would I have to define another
potential type (say Restraint and basically edit all code that deals with
the potential functions).
You can exclude pair types (not bond types) from pairwise interactions
via neigh_modify exclude. The delete_bonds command will also turn
off bonds of a certain type(s). But when they are off, the pairwise may
then turn on, depending on your special_bonds settings.
You could also create your own bond_ABC.cpp that computed each
bond type however you wanted.
Actually, I want to do the reverse. I have two bond types, a normal
harmonic one and a distance-based restraint (which is harmonic outside the
two bounds but 0 within the bounds). I want to exclude the harmonic bonds
(which are mapped onto real bonds) from pairwise interactions while
maintaining pair interactions between atoms that have a distance
restraint. Is this possible (or would it be difficult to add to the
code)? Looking at the neigh_modify command, basically I want "exclude"
based on bond type instead of atom type.
