Dear all,
I want to calculate mean squared displacements of atoms in CuZr using Einstein crystal model, which requires setting bond coefficients. From http://lammps.sandia.gov/threads/msg12670.html I learned that this can be done by setting the bond between each atom and its (equilibrium) lattice site via the so-called dummy atom method. But I could not find any information on how to set some dummy sites. Could someone help me? Thanks.
Eugene
Dear all,
I want to calculate mean squared displacements of atoms in CuZr using
Einstein crystal model, which requires setting bond coefficients. From
http://lammps.sandia.gov/threads/msg12670.html I learned that this can
be done by setting the bond between each atom and its (equilibrium)
lattice site via the so-called dummy atom method. But I could not find
any information on how to set some dummy sites. Could someone help me?
a dummy atom for your purpose would just be like a regular
atom that would be placed on the desired lattice position, but
then in LAMMPS it would not be included in the time integration.
LAMMPS is not directly expecting users to run this kind of
calculation, buy with a little experimentation and some creative
tricks it should be possible.
what you would need are the following steps:
- create a crystal of reference or dummy atom positions
assign them atom type 3 and assign them a group
- create a lattice each of Cu and Zr atom and assign them
atom types 1 and 2 respectively and assign them to a second group
- create a harmonic bond between each Cu or Zr atom and its
corresponding reference site and give it bond type 1 and 2, respectively
- use neigh_modify to exclude all pair interactions. you probably need
to assign some pair potential parameters, just us the lj/cut style and
set the coefficients * * 0.0 0.1
because of the bonds, you'll need to input the system from
a data file. if you can compute the lattice coordinate and write
them to an .xyz or .pdb file, you can build that in VMD using a
plugin that i wrote for it.
http://sites.google.com/site/akohlmey/software/topotools
see step1d of the associated tutorial about how to set arbitrary bonds:
https://sites.google.com/site/akohlmey/software/topotools/topotools-tutorial---part-1
let us know if you need more help.
axel.
Dear Axel (and others),
Thanks. I performed the first 3 steps you mentioned below and added line "atom_modify sort 0 0.0" which disables atom sorting. Since only harmonic bonding is specified, I assume pair interactions are excluded. Now a new problem arises.
I fixed the temperature of (Cu,Zr) group atoms to 500K, but at step 0, the temperature was 250K (according to thermo_style output). This was because the dummy atoms have 0 kinetic energy and temperature and thus decrease the temp of the whole system. If I used "thermo_modify temp ctempcz", where ctempcz represents a temperature calculation of (Cu,Zr) group, the temp at step 0 became 500 K, but the temp change towards to 1000 K as the system relaxes towards equilibrium. So it appears to me that the dummy atoms were included during time integration, which was confirmed by that the potential energy of (Cu,Zr) group was 1/2 that of the total system. Could you give me some hints how to exclude dummy atoms during time integration?
This Einstein crystal calculation is going to be a reference state for later free energy calculation using "compute ti" in LAMMPS. So I guess this should be solvable. Thanks.
Eugene
Quoting Axel Kohlmeyer <[email protected]>:
Dear Axel (and others)
I forgot to mention that during the calculation, the kinetic E of the dummy group is always 0 as step increases. This seems consistent with the requirement that dummy atoms be excluded. However, since Einstein model is harmonic, it's confusing that the potential energy of the dummy group (Ep_dm below) never converts into its kinetic energy (Ek_dm). A few lines from the log file is shown here for your reference. Some digits after the decimal point are removed for clarity and the subscript cz means (Cu,Zr) group.
Step Temp Temp_cz PotEng Ep_cz Ep_dm KinEng Ek_cz Ek_dm TotEng
0 500 500 0 0 0 9.2 9.2 0 9.2
100 12.8 12.8 10.1 5.0 5.0 0.2 0.2 0 10.3
...
9900 1075.0 1075.0 12.2 6.1 6.1 19.8 19.8 0 32.0
10000 868.3 868.3 16.1 8.0 8.0 16.0 16.0 0 32.2
Thanks. Eugene
Quoting [email protected]...:
Dear Axel (and others),
Thanks. I performed the first 3 steps you mentioned below and added
line "atom_modify sort 0 0.0" which disables atom sorting. Since only
harmonic bonding is specified, I assume pair interactions are
excluded. Now a new problem arises.I fixed the temperature of (Cu,Zr) group atoms to 500K, but at step 0,
the temperature was 250K (according to thermo_style output). This was
because the dummy atoms have 0 kinetic energy and temperature and thus
decrease the temp of the whole system. If I used "thermo_modify temp
ctempcz", where ctempcz represents a temperature calculation of
(Cu,Zr) group, the temp at step 0 became 500 K, but the temp change
towards to 1000 K as the system relaxes towards equilibrium. So it
appears to me that the dummy atoms were included during time
integration, which was confirmed by that the potential energy of
(Cu,Zr) group was 1/2 that of the total system. Could you give me some
hints how to exclude dummy atoms during time integration?
since i (and others) cannot read minds of people or computers,
please post your input to the list. otherwise it is near impossible
to say what you have done right and what is wrong.
axel.
Dear Axel,
Please see the attached in.lmp and geoin.lmp files. I used LAMMPS (15 Jan 2011). Thanks.
Eugene
Quoting Axel Kohlmeyer <[email protected]>:
in.lmp (1.19 KB)
geoin.lmp (27.5 KB)
Dear Axel (and others),
Thanks. I performed the first 3 steps you mentioned below and added
line "atom_modify sort 0 0.0" which disables atom sorting. Since only
harmonic bonding is specified, I assume pair interactions are
excluded. Now a new problem arises.I fixed the temperature of (Cu,Zr) group atoms to 500K, but at step 0,
the temperature was 250K (according to thermo_style output). This was
because the dummy atoms have 0 kinetic energy and temperature and thus
decrease the temp of the whole system. If I used "thermo_modify temp
ctempcz", where ctempcz represents a temperature calculation of
ok. that is reasonable. step, but you are missing some other issue:
you have to initialize your system to 1000K, since you are starting
all oscillators in the minimum. at equilibrium, you will have half of
that kinetic energy transferred into potential energy.
i also, don't understand why you are using an npt integrator.
you are defining your idealized lattice with your dummy atoms
and don't compute any interactions except the individual bonds.
so different particles don't really "see" each other.
(Cu,Zr) group, the temp at step 0 became 500 K, but the temp change
towards to 1000 K as the system relaxes towards equilibrium. So it
actually, there is not much equilibration going on. you can see that, if you
print out thermodynamic properties more often. the various oscillators
are happily wiggling around their equilibrium. since they all started in the
minimum, they are well correlated and only the coupling through the
thermostat algorithm will - very slowly - exchange energy between them
(indirectly via exchange with the "heat bath")/
appears to me that the dummy atoms were included during time
integration, which was confirmed by that the potential energy of
no they are not.
(Cu,Zr) group was 1/2 that of the total system. Could you give me some
that is how the energy is divided in bonds. each atom gets half.
hints how to exclude dummy atoms during time integration?
it was done correctly.
since i (and others) cannot read minds of people or computers,
please post your input to the list. otherwise it is near impossible
to say what you have done right and what is wrong.axel.
This Einstein crystal calculation is going to be a reference state for
later free energy calculation using "compute ti" in LAMMPS. So I guess
this should be solvable. Thanks.
another mystery to me is why you are computing the MSD.
the atoms cannot move around, they'll just oscillate around
their equilibrium position and not move anywhere.
axel.
One other option I forgot about in the code, which might
simplify your goal of tethering atoms to their initial site
with a harmonic spring: check out the fix spring/self
command. If this does what you want, you don't need
dummy atoms and to explicitly create bonds.
Steve
Dear Axel (and others)
After some thoughts, I did the following: (1) fix the (Cu,Zr) group using nve, which is more reasonable than npt for a harmonic oscillator. I used npt before just wanting to be consistent with the target system (Cu-Zr alloy using eam potential and npt, which corresponds to lambda=0 in "compute ti" method). (2) correspondingly, I only specified the initial velocity distribution (instead of fix temp) corresponding to 1000 K so that average E_kin and E_pot equals 500 K. This fixed my problem.
To your question why I want to calculate msd of Einstein crystal since it's only a harmonic oscillator: My goal is to get free energy of the Cu-Zr system using thermodynamical integration (using Einstein crystal as reference state and do "compute ti"). This reference state is preferably as close as possible to the real target system, which usually is determined by their msd. From msd of the target system, I can find the spring constant of the Einstein crystal by assuming that these two systems have the same msd. So I calculate the msd of the latter just for a confirmation.
Eugene
Quoting Axel Kohlmeyer <[email protected]>:
Dear Steve,
Yes. After a few try-and-error, I was able to do it, provided that I use "atomic" atom_style. So my related settings are:
atom_style atomic
fix fscu cu spring/self 1.172959
fix_modify fscu energy yes
fix fszr zr spring/self 1.125567
fix_modify fszr energy yes
But with this, I am confused about how to set the "compute ti" computation. Some of my related settings are:
variable uh equal lambda # related Ep is lambda*0.5*spring*(dr)^2
variable udh equal 1.0
variable ueam equal 1.0-lambda # related Ep is (1-lambda)*(Ep from eam potential)
variable udeam equal -1.0
compute cti_cz cz ti harmonic v_uh v_duh eam/fs v_ueam v_dueam
I knew that "harmonic" in the above "compute" line is not correct, and indeed I got an error message saying "ERROR: Variable for compute ti is invalid style". But I don't know what potential name I should use here as potential energy from "spring/self" is not defined(?). I also don't know if I need to declare variable "lambda" and explicitly specify its values like 0.1, 0.2,... A complete input file for "compute ti" is attached. Thanks.
Eugene
Quoting Steve Plimpton <[email protected]>:
in.ti (1.24 KB)
geoin.lmp-self (28.2 KB)
You are correct that you are not using the compute ti command
correctly - I would read the doc page carefully. It is setup to
use with pair styles and Kspace. Fix spring/self is neither of those,
as it is not a potential and it is not controllable by lamda. You may
want to talk to Sai (CCd) and see if he has ideas about computing
the free energy of your system. He is the expert here on that.
Steve