Hello everybody,
I am trying to set rigid molecules in my simulation. I tried using fix shake and fix rigid. In case of fix shake it asked me for bond potentials which I dont have for the molecule. Whereas, in case of fix rigid the energies are highly positive. It would be highly appreciable if anyone could suggest some clue to the solution. Thanks in advance.
Kind regards
Bhaskar
Hello everybody,
I am trying to set rigid molecules in my simulation. I tried using fix shake
and fix rigid. In case of fix shake it asked me for bond potentials which I
dont have for the molecule. Whereas, in case of fix rigid the energies are
you have a molecule without bonds? how is this possible??
highly positive. It would be highly appreciable if anyone could suggest some
clue to the solution. Thanks in advance.
the only clue that i can give at this point is that _you_ have to
give _us_ (i.e. people on the mailing list) more of a clue about
what it is that you are doing. there are a ton of things that can
lead to high energies, most of which are mistakes on your side.
yet without knowing whether you are trying to simulate a family
of elephants, some out-of-space compound or simply plain water,
there is little help to be offered.
axel.
Hi Axel.
Here is the details of what I am trying to do. I am trying to simulate a mixture of methane and CO2 at normal temperature and pressure. The CH4 is an united atom single bead in the model I am using. On the other hand CO2 is a linear tri atomic molecule with it bonds and angle fixed. Basically it is a rigid molecule. In the CO2 model I do not have the force constants for bond vibration as well as bending. The initial configuration was taken from a short TOWHEE run.
When I use shake to constrain the bonds in CO2, LAMMPS ask me for the bond potential.
I also tried with fix rigid. In this case the simulation runs but the energy is highly positive. Below I have pasted the input script. If you could kindly go through it and comment, it will be really helpful.
units real
atom_style full
boundary p p p
pair_style lj/cut/coul/cut 13.3 13.3
pair_modify shift yes mix arithmetic
#kspace_style pppm 1.0e-4
dielectric 1.0
read_data co2_ch4.data
#Groups for the molecules
group CO2 type 2 3
group Methane type 1
timestep 1.0 #femto seconds
velocity all create 293.15 2309483 mom yes rot yes dist gaussian
fix NVE Methane nve
fix RIGID_CO2 CO2 rigid molecule
fix TEMPRESCALE all temp/rescale 1 293.15 293.15 0.0 1.0
compute TEMP_all all temp
thermo 1000
thermo_style custom step c_TEMP_all press vol ke pe etotal ecoul
thermo_modify flush yes
#Initial equilibration steps
dump ATOM_COOR all custom 1000 dump.atom id mol type x y z vx vy vz
run 100000 #Equilibration steps
kind regards
Bhaskar
Hi Axel.
Here is the details of what I am trying to do. I am trying to simulate a
mixture of methane and CO2 at normal temperature and pressure. The CH4 is an
united atom single bead in the model I am using. On the other hand CO2 is a
linear tri atomic molecule with it bonds and angle fixed. Basically it is a
rigid molecule. In the CO2 model I do not have the force constants for bond
vibration as well as bending. The initial configuration was taken from a
short TOWHEE run.
well, the value of the force constant has no meaning when using shake,
so this would not be a problem, but - as far as i remember - the shake
implementation in LAMMPS cannot handle linear molecules, so it is
pretty much pointless to try using it in any case. you should search
the mailing list archives, since this topic has been discussed extensively
a number of times in the past.
When I use shake to constrain the bonds in CO2, LAMMPS ask me for the bond
potential.
yes. but all it uses is the bond length.
I also tried with fix rigid. In this case the simulation runs but the energy
is highly positive. Below I have pasted the input script. If you could
kindly go through it and comment, it will be really helpful.
the input script alone is usually not sufficient. one can only run a proper
test and evaluate the impact of changes to various settings by having a
_complete_ input (i.e. including the data file).
however, there are two items in your input script that are worrying me:
you using a 1.0 fs time step and you are using temp/rescale.
fix rigid integrates equations of motion differently than when using
fix nve with fix shake and as a consequence you need to have a
shorter time step usually. and fix temp/rescale is in general a
horrible way to "control" temperature, since it is the same as
kicking the system all the time. there are fix rigid/nvt and fix nvt
or you can use fix langevin instead of fix temp/rescale.
axel.