enthalpy and equation of state

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 --------------------- 


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 


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



print "creo ya esta =)"

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?

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, …

Oscar Guerrero

The loop is easy to implement using a loop type variable in lammps.
The variable link show one example: