I am a new lammps user (rpm version 4.11.2, Jul 25 2016). I have a
problem with my first script. I want to calculate lattice parameter
as a function of temperature (1 step) and then the stress generated by rapid cooling (2 step). For this I have to stabilize system at high temperature - especially pressure (1 step). And here (1 step) I have a problem with the pressure fluctuation at high temperature that are very much higher than the temperature fluctuation. I attach my script (for 1 step) for a constant temperature (600 K) and pressure (1 bar). Attached script is for Al because result can be compared with the experimental data from the literature. My goal is a intermetalic compound that does not have such data.
I know that two fixes are not physical but fix2 (when active) significantly reduces pressure fluctuations. Can I do anything more or increase the system is the only way to improve? I will also be grateful for other corrections, especially information how to enter internal fix 1 data (temperature, pressure, volume) into output file.
Greetings
Mirek
-----------------------
units metal
boundary p p p
atom_style atomic
pair_style meam
lattice fcc 4.078
region box block 0 16 0 16 0 16
create_box 1 box
lattice fcc 4.078 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
create_atoms 1 box
mass 1 26.98
pair_style eam/alloy
pair_coeff * * Al99.eam.alloy Al
neighbor 2.0 bin
neigh_modify delay 0 every 1
compute eng all pe/atom
compute eatoms all reduce sum c_eng
reset_timestep 0
timestep 0.0001
velocity all create 1000 12345678 mom yes rot no
fix 1 all npt temp 600 600 1.0 iso 1.0 1.0 1.20 drag 0.5
# fix 2 all press/berendsen iso 1.0 1.0 10000.0
compute 2 all temp
thermo 500
thermo_style custom step temp press vol lx xlat
minimize 1e-30 1e-30 1000000 1000000
run 1000000
unfix 1
# unfix 2
variable natoms equal "count(all)"
variable teng equal "c_eatoms"
variable length equal "lx/10"
variable ecoh equal "v_teng/v_natoms"
print "Total energy (eV) = ${teng};"
print "Number of atoms = ${natoms};"
print ${length}
print "Lattice constant (Angstoms) = ${length};"
print "Cohesive energy (eV) = ${ecoh};"
print "All done!".
Your Tdamp and Pdamp parameters are quite high. Values equal to 100 timesteps for Tdamp and 1000 timesteps for Pdamp should work better. In order to lower the initial pressure add fix box/relax before minimization.
You can output the data using fix print with "file" keyword.
I am not familiar with this software. But, I would like to make one remark. NpT is useful as long as you do not believe to get the real pressure using the virial in MD or MC. There are several sources of error particularly for nontrivial molecular structures including partial charges on the atoms. For checking the temperature dependence of a lattice size under constant 'pressure' it might be o.k. if you do not need the exact value of the pressure.
There is at least one issue with your script.
Line: velocity all create 1000 12345678 mom yes rot no ; says heat to 1000K
Following line: fix 1 all npt temp 600 600 1.0 iso 1.0 1.0 1.20 drag 0.5 ; says that the initial temperature is 600 K, so from one line to the next one the program is trying to adjust for a difference of 400K.
My suggestion to reduce fluctuations is:
1) Perform a short NVE run after the velocity create command, check pressure of your sample.
2) Perform NPT at constant temperature (1000K) beginning with the final pressure of the NVE run to the target pressure (1 bar).
3) Perform NPT at constant pressure (1 bar) from 1000K to 600K.
4) Perform NPT at constant pressure and constant temperature.
Also, make sure the iso constrain is valid for your system under study.
You should see less fluctuations now.
Your Tdamp and Pdamp parameters are quite high. Values equal to 100 timesteps for Tdamp and 1000 timesteps for Pdamp should work better. In order to lower the initial pressure add fix box/relax before minimization.
You can output the data using fix print with "file" keyword.
change your sript:
velocity all create 1000 12345678 mom yes rot no
fix 1 all npt temp 600 600 1.0 iso 1.0 1.0 1.20 drag 0.5
# fix 2 all press/berendsen iso 1.0 1.0 10000.0
compute 2 all temp
thermo 500
thermo_style custom step temp press vol lx xlat
minimize 1e-30 1e-30 1000000 1000000
run 1000000
as:
fix 100 all box/relax iso 1 vmax 0.1
minimize 1e-30 1e-30 1000000 1000000
velocity all create 600 12345678 mom yes rot no
fix 1 all npt temp 600 600 0.01 iso 1.0 1.0 0.1
compute 2 all temp
thermo 500
thermo_style custom step temp press vol lx xlat
run 1000000
just try it, hope it can help you
It is important to reduce the forces by static equilibration of the system, before to start the MD simulation (by command minimize). Otherwise the initial velocities are much higher than those corresponding to the target temperature and the MD equilibration stage will be very long. You also may eliminate the rotation (by command velocity all create 600 12345678 mom yes rot yes) of the system.
Dear colleagues,
I answer only now because your proposition has been tested by me. Relaxation proposed by Xin He works very well. Explanations given by Carlos Ruestes also helped me to improve the program. I tested a damping parameters given by Michal Kanski. They are the same as in the LAMMPS manual. However, some larger damping values give smaller oscillations in my case. Letters sent by Siegfrid Fritzsche and the Vioriel Chihaia were also very helpful.
Thank you very much for your help in solving my problem.
Greetings,
M. Wrobel
Dear Michal Kanski
In your answer you mentioned: 100 timesteps for Tdamp and 1000 timesteps for Pdamp
suppose my time step is 0.002
then whether Tdamp is 100*0.002 or just 100. and for pdamp is it 1000*0.002 or just 1000 irrespective of time step. can you give full details?
- Badges
- Science topic