I am trying to minimize structure which is a hairpin-like structure (input structure inluded as system.data). There is no solvent in the system.
Energy minimization gives me -nan for E_mol and pressure is also nan. I have included the input file that I am using here:
# -- initialise --
units real
atom_style full
boundary p p p
timestep 10
# -- read input --
read_data "system.data"
# -- force field --
include "ff.settings"
mass 1 14.0270
set type 1 charge 0
# -- group --
group pol type 1
# -- simulation protocol --
# EM
timestep 2
thermo 1
minimize 0.0 0.0 1000 10000
write_data postem.data
The ouput I get is:
Step Temp E_pair E_mol TotEng Press
0 0 229.63578 -nan -nan nan
1 0 229.63578 -nan -nan nan
I though energy minimisation should take care of the unphysical structure I have by automatically making the angles and bond lengths to the correct values. Could someone tell me what am I missing ?
Edit: added force field file ff.settings
Thanks and regards.
Depending on the analytical form of your molecular interactions, some geometries can lead to mathematically undefined behavior. Many angular energy forms in cos(\theta), either proper, dihedral or improper, lead to division by sin(\theta) or other oddities when computing the forces. Obviously special values like \theta=\pi can result in divergences as mentioned by @akohlmey.
Since the energy term resulting in NaN is the molecular interaction, I also suspect your issue stems from something like that. Adding some noise or using values just aside the problematic values often correct these problems.
I do have cosine in angle potential because I want to use cosine squared angle potential.
So are you suggesting that the initial geometry is problematic? If I try to displace atoms by using displace_atoms and displace them by 0.5 or 1.0, the output does not change
The only term that could cause that is the angle. By doing the minimisation with angle_style zero and then reverting to the intended angle style, the simulation runs just fine.
polymer.in
# -- initialise --
units real
atom_style full
boundary p p p
# -- read input --
read_data "system.data"
# Forcefield parameters
pair_style lj/cut 10.0
# eps sig
pair_coeff 1 1 0.239006 4.00
bond_style harmonic
bond_coeff 1 394.3595 1.530 # P-P
# Just for the minimisation
angle_style zero
angle_coeff 1 11.9503
special_bonds lj 0 0 0.5
set type 1 charge 0
# -- group --
group pol type 1
# -- simulation protocol --
dump 1 all atom 10 polymer.dump
thermo 1
minimize 0.0 0.0 1000 10000
write_data postem.data
# Back to the full potential.
angle_style cosine/squared
angle_coeff 1 11.9503 111.00 # P-P-P
# Short annealing.
timestep 2
velocity all create 300 366929 dist gaussian rot yes
fix 1 all nvt temp 300 300 $(100*dt)
run 100000
Yes indeed, it works. However, the structure I get has some wierd helix in the bend of the polymer which is not exactly what I want.
However, if the angle potential being cosine is a problem, then making the angle potential harmonic should also be fine right? However, if I make it harmonic, I still get -nan for E_mol.
The problem is the initial structure. You should change the angle potential after the energy minimisation with style zero.
Regarding the desired morphology, it is your job to model correctly the interactions in your system. LAMMPS only integrates the equations of motions for you and the system follows the physics defined by the fixes and force field parameters. There is no much to add, but a nice animation for the other users
