Skip to content

Commit

Permalink
Merge pull request #13 from ICAMS/update/liquid
Browse files Browse the repository at this point in the history 
Update/liquid
  • Loading branch information
@srmnitc
srmnitc authored Nov 3, 2021
2 parents ffe71b1 + 9f1dba4 commit edee67a
Show file tree
Hide file tree
Showing 9 changed files with 338 additions and 417 deletions.
2 changes: 1 addition & 1 deletion .bumpversion.cfg
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
[bumpversion]
current_version = 0.5.8
current_version = 0.6.0
commit = True
tag = True

Expand Down
14 changes: 12 additions & 2 deletions README.md
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,16 @@
# calphy

Python library and command line tool for calculation of free energies.
`calphy`(pronounced _cal-phee_) is a Python library and command line tool for calculation of free energies. It uses [LAMMPS](https://www.lammps.org/) as the molecular dynamics driver to enable calculation of free energies using thermodynamic integration in a completely automated and efficient manner. The various methods that calphy can perform are:

- $F(V_i,T_i)$ and $G(P_i,T_i)$ for both solid and liquid phases at the given thermodynamic conditions using [non-equilibrium Hamiltonian interpolation](https://linkinghub.elsevier.com/retrieve/pii/S0927025615007089).
- $F(T_i \to T_f)$ and $G(T_i \to T_f)$, temperature dependence of Gibbs and Helmholtz free energies using the [reversible scaling](https://link.aps.org/doi/10.1103/PhysRevLett.83.3973) approach.
- Calculation of solid-solid or solid-liquid phase transition temperatures.
- Calculation of coexistence lines using [dynamic Clausius-Clapeyron integration](http://aip.scitation.org/doi/10.1063/1.1420486).
- Calculation of specific heat $c_P(T)$ as a function of temperature.
- Calculation of $F(x, T)$ and $G(x, T)$ for multicomponent systems using [alchemical transformations](https://journals.aps.org/prmaterials/abstract/10.1103/PhysRevMaterials.5.103801).
- [Upsampling calculations](https://journals.aps.org/prmaterials/abstract/10.1103/PhysRevMaterials.5.103801).

Calphy works with all [interatomic potentials implemented in LAMMPS](https://docs.lammps.org/pairs.html) and can obtain reliable results with low error bars in as low as 50 ps of simulation time with less than 3000 atoms. More information about the methods in calphy can be found in the [associated publication](https://journals.aps.org/prmaterials/abstract/10.1103/PhysRevMaterials.5.103801).

## Installation, examples and documentation

Expand All @@ -13,7 +23,7 @@ Python library and command line tool for calculation of free energies.
## Citing calphy

If you find calphy useful, please consider citing:
Menon, Sarath, Yury Lysogorskiy, Jutta Rogal, and Ralf Drautz.
Menon, Sarath, Yury Lysogorskiy, Jutta Rogal, and Ralf Drautz.
[“Automated Free Energy Calculation from Atomistic Simulations.”](https://journals.aps.org/prmaterials/abstract/10.1103/PhysRevMaterials.5.103801)
Physical Review Materials 5(10), 2021
DOI: 10.1103/PhysRevMaterials.5.103801
Expand Down
12 changes: 6 additions & 6 deletions calphy/integrators.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -229,12 +229,12 @@ def integrate_path(fwdfilename, bkdfilename, nelements=1, concentration=[1,], us
# - Liquid : UFM system energy is not scaled with fix:adapt (this is actually wrong, but works for us)
# - Liquid : So we only unscale the system of interest in liquid
# - Alchemy : Both system energies are scaled with fix:adapt, so we need to unscale both
for i in range(len(fdui)):
if flambda[i] !=0:
fdui[i] = fdui[i]/flambda[i]
for i in range(len(bdui)):
if blambda[i] !=0:
bdui[i] = bdui[i]/blambda[i]
#for i in range(len(fdui)):
# if flambda[i] !=0:
# fdui[i] = fdui[i]/flambda[i]
#for i in range(len(bdui)):
# if blambda[i] !=0:
# bdui[i] = bdui[i]/blambda[i]

if alchemy:
"""
Expand Down
274 changes: 104 additions & 170 deletions calphy/liquid.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -228,76 +228,120 @@ def run_integration(self, iteration=1):
lmp = ph.read_dump(lmp, conf, species=self.options["nelements"])

#set hybrid ufm and normal potential
lmp = ph.set_hybrid_potential(lmp, self.options, self.eps)
#lmp = ph.set_hybrid_potential(lmp, self.options, self.eps)
lmp = ph.set_potential(lmp, self.options)

#remap the box to get the correct pressure
lmp = ph.remap_box(lmp, self.lx, self.ly, self.lz)

#apply the necessary thermostat
lmp.command("fix f1 all nve")
lmp.command("fix f2 all langevin %f %f %f %d"%(self.t, self.t, self.options["md"]["tdamp"],
np.random.randint(0, 10000)))

lmp.command("fix f1 all nve")
lmp.command("fix f2 all langevin %f %f %f %d zero yes"%(self.t, self.t, self.options["md"]["tdamp"],
np.random.randint(0, 10000)))
lmp.command("run %d"%self.options["md"]["te"])

# Compute the potential energy of each pair style.
lmp.command("compute c1 all pair %s"%self.options["md"]["pair_style"])
lmp.command("compute c2 all pair ufm")
lmp.command("unfix f1")
lmp.command("unfix f2")

# Output variables.
lmp.command("variable step equal step")
lmp.command("variable dU1 equal c_c1/atoms") # Driving-force obtained from NEHI procedure.
lmp.command("variable dU2 equal c_c2/atoms")
#---------------------------------------------------------------
# FWD cycle
#---------------------------------------------------------------

#force thermo to evaluate variables
lmp.command("thermo_style custom step v_dU1 v_dU2")
lmp.command("thermo 1000")
lmp.command("variable flambda equal ramp(${li},${lf})")
lmp.command("variable blambda equal 1.0-v_flambda")

#switching completely to potential of interest
lmp.command("variable zero equal 0")
lmp.command("fix f0 all adapt 0 pair ufm scale * * v_zero")
lmp.command("run 0")
lmp.command("unfix f0")
lmp.command("pair_style hybrid/scaled v_flambda %s v_blambda ufm 7.5"%self.options["md"]["pair_style"])

#Equilibrate system
lmp.command("run %d"%self.options["md"]["te"])
pc = self.options["md"]["pair_coeff"]
pcraw = pc.split()
pcnew = " ".join([*pcraw[:2], *[self.options["md"]["pair_style"],], *pcraw[2:]])

#print header
lmp.command("print \"${dU1} ${dU2} ${li}\" file forward_%d.dat"%iteration)

#set up scaling variables
lmp.command("variable lambda_p1 equal ramp(${li},${lf})")
lmp.command("variable lambda_p2 equal ramp(${lf},${li})")

#Forward switching run
lmp.command("fix f3 all adapt 1 pair %s scale * * v_lambda_p1"%self.options["md"]["pair_style"])
lmp.command("fix f4 all adapt 1 pair ufm scale * * v_lambda_p2")
lmp.command("fix f5 all print 1 \"${dU1} ${dU2} ${lambda_p1}\" screen no append forward_%d.dat"%iteration)
lmp.command("run %d"%self.options["md"]["ts"])

#unfix things
lmp.command("unfix f3")
lmp.command("unfix f4")
lmp.command("unfix f5")

#Equilibriate at UFM potential
lmp.command("run %d"%self.options["md"]["te"])

#print file header
lmp.command("print \"${dU1} ${dU2} ${lf}\" file backward_%d.dat"%iteration)

#set up scaling variables
lmp.command("variable lambda_p1 equal ramp(${lf},${li})")
lmp.command("variable lambda_p2 equal ramp(${li},${lf})")

#Reverse switching run
lmp.command("fix f3 all adapt 1 pair %s scale * * v_lambda_p1"%self.options["md"]["pair_style"])
lmp.command("fix f4 all adapt 1 pair ufm scale * * v_lambda_p2")
lmp.command("fix f5 all print 1 \"${dU1} ${dU2} ${lambda_p1}\" screen no append backward_%d.dat"%iteration)
lmp.command("run %d"%self.options["md"]["ts"])

#unfix things
lmp.command("unfix f3")
lmp.command("unfix f4")
lmp.command("unfix f5")
lmp.command("pair_coeff %s"%pcnew)
lmp.command("pair_coeff * * ufm %f 1.5"%self.eps)

lmp.command("compute c1 all pair %s"%self.options["md"]["pair_style"])
lmp.command("compute c2 all pair ufm")

lmp.command("variable step equal step")
lmp.command("variable dU1 equal c_c1/atoms")
lmp.command("variable dU2 equal c_c2/atoms")

lmp.command("thermo_style custom step v_dU1 v_dU2")
lmp.command("thermo 1000")


lmp.command("velocity all create %f %d mom yes rot yes dist gaussian"%(self.t, np.random.randint(0, 10000)))

lmp.command("fix f1 all nve")
lmp.command("fix f2 all langevin %f %f %f %d zero yes"%(self.t, self.t, self.options["md"]["tdamp"],
np.random.randint(0, 10000)))
lmp.command("compute Tcm all temp/com")
lmp.command("fix_modify f2 temp Tcm")

lmp.command("fix f3 all print 1 \"${dU1} ${dU2} ${flambda}\" screen no file forward_%d.dat"%iteration)
lmp.command("run %d"%self.options["md"]["ts"])

lmp.command("unfix f1")
lmp.command("unfix f2")
lmp.command("unfix f3")
lmp.command("uncompute c1")
lmp.command("uncompute c2")

#---------------------------------------------------------------
# EQBRM cycle
#---------------------------------------------------------------

lmp.command("pair_style ufm 7.5")
lmp.command("pair_coeff * * %f 1.5"%self.eps)

lmp.command("thermo_style custom step pe")
lmp.command("thermo 1000")

lmp.command("fix f1 all nve")
lmp.command("fix f2 all langevin %f %f %f %d zero yes"%(self.t, self.t, self.options["md"]["tdamp"],
np.random.randint(0, 10000)))
lmp.command("fix_modify f2 temp Tcm")

lmp.command("run %d"%self.options["md"]["te"])

lmp.command("unfix f1")
lmp.command("unfix f2")

#---------------------------------------------------------------
# BKD cycle
#---------------------------------------------------------------

lmp.command("variable flambda equal ramp(${lf},${li})")
lmp.command("variable blambda equal 1.0-v_flambda")

lmp.command("pair_style hybrid/scaled v_flambda %s v_blambda ufm 7.5"%self.options["md"]["pair_style"])

lmp.command("pair_coeff %s"%pcnew)
lmp.command("pair_coeff * * ufm %f 1.5"%self.eps)

lmp.command("compute c1 all pair %s"%self.options["md"]["pair_style"])
lmp.command("compute c2 all pair ufm")

lmp.command("variable step equal step")
lmp.command("variable dU1 equal c_c1/atoms")
lmp.command("variable dU2 equal c_c2/atoms")

lmp.command("thermo_style custom step v_dU1 v_dU2")
lmp.command("thermo 1000")

lmp.command("fix f1 all nve")
lmp.command("fix f2 all langevin %f %f %f %d zero yes"%(self.t, self.t, self.options["md"]["tdamp"],
np.random.randint(0, 10000)))
lmp.command("fix_modify f2 temp Tcm")

lmp.command("fix f3 all print 1 \"${dU1} ${dU2} ${flambda}\" screen no file backward_%d.dat"%iteration)
lmp.command("run %d"%self.options["md"]["ts"])

lmp.command("unfix f1")
lmp.command("unfix f2")
lmp.command("unfix f3")
lmp.command("uncompute c1")
lmp.command("uncompute c2")

#close object
lmp.close()
Expand Down Expand Up @@ -347,113 +391,3 @@ def thermodynamic_integration(self):
self.fe = self.fideal + self.fref - self.w + self.pv



def reversible_scaling(self, iteration=1):
"""
Perform reversible scaling calculation in NPT
Parameters
----------
iteration : int, optional
iteration of the calculation. Default 1
Returns
-------
None
"""
t0 = self.t
tf = self.tend
li = 1
lf = t0/tf
pi = self.p
pf = lf*pi

#create lammps object
lmp = ph.create_object(self.cores, self.simfolder, self.options["md"]["timestep"])

lmp.command("echo log")
lmp.command("variable li equal %f"%li)
lmp.command("variable lf equal %f"%lf)

#read in conf file
conf = os.path.join(self.simfolder, "conf.dump")
lmp = ph.read_dump(lmp, conf, species=self.options["nelements"])

#set up potential
lmp = ph.set_potential(lmp, self.options)

#remap the box to get the correct pressure
lmp = ph.remap_box(lmp, self.lx, self.ly, self.lz)

#set thermostat and run equilibrium
lmp.command("fix f1 all nph iso %f %f %f"%(self.p, self.p, self.options["md"]["pdamp"]))
lmp.command("fix f2 all langevin %f %f %f %d zero yes"%(self.t, self.t, self.options["md"]["tdamp"], np.random.randint(0, 10000)))
lmp.command("run %d"%self.options["md"]["te"])
lmp.command("unfix f1")
lmp.command("unfix f2")

#now fix com
lmp.command("variable xcm equal xcm(all,x)")
lmp.command("variable ycm equal xcm(all,y)")
lmp.command("variable zcm equal xcm(all,z)")

lmp.command("fix f1 all nph iso %f %f %f fixedpoint ${xcm} ${ycm} ${zcm}"%(self.p, self.p, self.options["md"]["pdamp"]))
lmp.command("fix f2 all langevin %f %f %f %d zero yes"%(t0, t0, self.options["md"]["tdamp"], np.random.randint(0, 10000)))

#compute com and modify fix
lmp.command("compute tcm all temp/com")
lmp.command("fix_modify f1 temp tcm")
lmp.command("fix_modify f2 temp tcm")

lmp.command("variable step equal step")
lmp.command("variable dU equal c_thermo_pe/atoms")
lmp.command("thermo_style custom step pe c_tcm press vol")
lmp.command("thermo 10000")

#create velocity and equilibriate
lmp.command("velocity all create %f %d mom yes rot yes dist gaussian"%(t0, np.random.randint(0, 10000)))
lmp.command("run %d"%self.options["md"]["te"])

#unfix nph
lmp.command("unfix f1")

#define ramp
lmp.command("variable lambda equal ramp(${li},${lf})")

#start scaling over switching time
lmp.command("fix f1 all nph iso %f %f %f fixedpoint ${xcm} ${ycm} ${zcm}"%(pi,
pf, self.options["md"]["pdamp"]))
lmp.command("fix_modify f1 temp tcm")
lmp.command("fix f3 all adapt 1 pair %s scale * * v_lambda"%self.options["md"]["pair_style"])
lmp.command("fix f4 all print 1 \"${dU} $(press) $(vol) ${lambda}\" screen no file forward_%d.dat"%iteration)
lmp.command("run %d"%self.options["md"]["ts"])

#unfix
lmp.command("unfix f3")
lmp.command("unfix f4")
lmp.command("unfix f1")

#equilibriate scaled hamiltonian
lmp.command("fix f1 all nph iso %f %f %f fixedpoint ${xcm} ${ycm} ${zcm}"%(pf,
pf, self.options["md"]["pdamp"]))
lmp.command("fix_modify f1 temp tcm")
lmp.command("run %d"%self.options["md"]["te"])
lmp.command("unfix f1")

#reverse lambda ramp
lmp.command("variable lambda equal ramp(${lf},${li})")

#apply fix and perform switching
lmp.command("fix f1 all nph iso %f %f %f fixedpoint ${xcm} ${ycm} ${zcm}"%(pf,
pi, self.options["md"]["pdamp"]))
lmp.command("fix_modify f1 temp tcm")
lmp.command("fix f3 all adapt 1 pair %s scale * * v_lambda"%self.options["md"]["pair_style"])
lmp.command("fix f4 all print 1 \"${dU} $(press) $(vol) ${lambda}\" screen no file backward_%d.dat"%iteration)
lmp.command("run %d"%self.options["md"]["ts"])
lmp.command("unfix f3")
lmp.command("unfix f4")
lmp.command("unfix f1")

#close the object
lmp.close()

Loading

0 comments on commit edee67a

Please sign in to comment.