Hi lammps users,
im trying to calculate the equation of state (EOS) (T=0) for different structures eg. BCC, FCC, HCP. Please find the script below. I'm using fix deform and storing Etotal/atom and VOL/Vo.
The plot of (E vs V) I believe make sence because the minimum of energy is equal to the Ecohesive. Foe example lets use Cu EAM Potential Mishin Cu EAM1 PRB(2001)63:224106
I'm also trying to find Enthalpy vs Pressure for different phases. I'm trying to Compare FCC enthalpy vc HCP enthalpy at T=0 .
Unfourtounately my script is wrong. In order to compare two data sets, the data sets should be of same lenght and the values space evenly i.e Equal values of Pressure P correspong to different values of Enthalpy "H"
Example ENTHALPHY VS P (BCC) P=10 , H=1 ENTHALPHY VS P (HCP) P=10 , H=23 (values are spaces evenly in pressure)
Can someone please tell me how to find Enthalpy vs Pressure (T=0) for different Phases using Lammps ?
Any information will be grealy appreciate
Thanks alot for your help and advise
Oscar Guerrero
# ---------- Initialize Simulation ---------------------
clear
units metal
dimension 3
boundary p p p
atom_style atomic
atom_modify map array
# ---------- define variables -------------------------------
timestep 0.005
# ---------- Create Atoms ------------------------------------
lattice bcc 3.615 origin 0.0 0.0 0.0 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
# lattice fcc 3.615 origin 0.0 0.0 0.0 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
# lattice hcp 3.615 origin 0.0 0.0 0.0 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
region box block 0 5 0 5 0 5 units lattice
create_box 1 box
create_atoms 1 box
# ---------- Define Interatomic Potential ---------------------
pair_style eam/alloy
pair_coeff * * Cu_mishin1.eam.alloy Cu
# ---------- Define Settings ----------------------------------
compute eng all pe/atom
compute kiene all ke/atom
compute pre all stress/atom
compute eatoms all reduce sum c_eng
# ---------- Print Atoms ----------------------------------
reset_timestep 0
dump 1 all cfg 2 lattice.dump.*.cfg id type xs ys zs
run 1
undump 1
# ---------- FIND MINIMUM ENERGY LEVEL -------------------------
reset_timestep 0
fix 1 all box/relax iso 0.0 vmax 0.001
thermo 0
thermo_style custom step pe lx ly lz vol press pxx pyy pzz c_eatoms
min_style cg
minimize 1e-25 1e-25 5000 10000
unfix 1
#--------------STORE TEMPORARY VARIABLES ----------------------------
# Store final volume/atom for further calculations
variable tmp equal "vol/count(all)"
variable V0 equal ${tmp}
print "Volume per atom, V0: ${V0}"
# Store final length for further calculations
variable tmp equal "lx"
variable L0 equal ${tmp}
print "Initial Length, L0: ${L0}"
# ---------- DEFORMATION STEPS ARE HERE ----------------------------------
reset_timestep 0
# RUN SIMULATION AT 0K
fix 1 all deform 1 x erate -0.10 remap x
# ---------- FIND EQUATION OF STATE ------------------------------------
thermo 0
thermo_style custom step etotal enthalpy lx ly lz vol press pxx pyy pzz temp
variable natoms equal "count(all)" # NUMBER ATOMS
variable H equal "enthalpy/v_natoms" # ENTHALPY
variable strain equal "-(lx - v_L0)/v_L0" # STRAIN
variable shear equal "0.5*(pxx - 0.5*(pyy + pzz))/10000" # SHEAR CALCULATION
variable out1 equal "(vol/v_natoms)/v_V0" # VOLUME PER ATOM
variable out2 equal "etotal/v_natoms" # TOTAL ENERGY PER/ATOM
variable out3 equal "pxx/10000" # PRESSURE PXX (GPA UNITS)
variable out4 equal "pyy/10000" # PRESSURE PYY (GPA UNITS)
variable out5 equal "pzz/10000" # PRESSURE PZZ (GPA UNITS)
variable out6 equal "press/10000" # TENSOR PRESSURE (GPA UNITS)
variable out7 equal "temp" # TEMPERATURE
variable out8 equal "lx" # lx
variable out9 equal "ly" # ly
variable out10 equal "lz" # lz
fix def1 all print 2 "${strain} ${shear}" file shear.dat screen no
fix def2 all print 2 "${out1} ${out2} ${out3} ${out4} ${out5} ${out6} ${out7} ${out8} ${out9} ${out10} ${H} " file eos.dat screen no
fix def3 all print 2 "${out6} ${H} " file enthalpy.dat screen no
run 1200
unfix 1
unfix def1
unfix def2
unfix def3
# SIMULATION DONE
clear
print "creo ya esta =)"
Oscar,
If I understand correctly your problem with the H vs P plot arises
because the variable that you modify with fix deform is the cell shape
(volume), and thus you get different pressure values when considering
different phases. Right?
Have you tried instead to place the fix box/relax combined with
minimization in a loop where you set the pressure within the fix
box/relax and then extract your other variables from the thermo info?
Carlos
Hi Carlos,
Your right, The problem arise when i try to compare (H vs P) for different phases. The output data for the two phases are not evenly spaced, thus its impossible to compare the solid-solid transition I came with another approach that might work: How about Using xmgrace and then do interpolation/spline (cubic splne) to the data sets. This way i can define the low/hight limit and the length, Using this mehod then the sets will be of same size and then i can compare Hfcc - Hhcp and find the pressure § in which the phase transition occours i.e Hfcc - Hhcp = 0
I Read your comment and your make me realice that when using fix deform i’m comparing the two phases, but I created two different lattice and thus the value of pressure are not going to be equal (duh!) . Please Let me think about the FOR/Loop feature (I follow you idear) and hopefuly i may be able to implent an script, …
Thanks
Oscar Guerrero
Oscar
The loop is easy to implement using a loop type variable in lammps.
The variable link show one example:
http://lammps.sandia.gov/doc/variable.html
Carlos