diff --git a/.bumpversion.cfg b/.bumpversion.cfg index 8f60cdd..5f5fe49 100644 --- a/.bumpversion.cfg +++ b/.bumpversion.cfg @@ -1,5 +1,5 @@ [bumpversion] -current_version = 0.5.8 +current_version = 0.6.0 commit = True tag = True diff --git a/README.md b/README.md index 86e02d9..48caed0 100644 --- a/README.md +++ b/README.md @@ -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 @@ -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 diff --git a/calphy/integrators.py b/calphy/integrators.py index 66782b9..7982d6a 100644 --- a/calphy/integrators.py +++ b/calphy/integrators.py @@ -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: """ diff --git a/calphy/liquid.py b/calphy/liquid.py index eed1825..9056b1a 100644 --- a/calphy/liquid.py +++ b/calphy/liquid.py @@ -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() @@ -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() - diff --git a/calphy/phase.py b/calphy/phase.py index 8a348e7..ca0cea0 100644 --- a/calphy/phase.py +++ b/calphy/phase.py @@ -249,6 +249,148 @@ def submit_report(self): elif self.calc["mode"] == "mts": self.logger.info("- 10.1063/1.1420486") + + def reversible_scaling(self, iteration=1, solid=False): + """ + 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 %s %f %f %f"%(self.iso, 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 %s %f %f %f fixedpoint ${xcm} ${ycm} ${zcm}"%(self.iso, 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") + + + lmp.command("variable flambda equal ramp(${li},${lf})") + lmp.command("variable blambda equal ramp(${lf},${li})") + lmp.command("variable fscale equal v_flambda-1.0") + lmp.command("variable bscale equal v_blambda-1.0") + lmp.command("variable one equal 1.0") + + #set up potential + pc = self.options["md"]["pair_coeff"] + pcraw = pc.split() + pcnew1 = " ".join([*pcraw[:2], *[self.options["md"]["pair_style"],], "1", *pcraw[2:]]) + pcnew2 = " ".join([*pcraw[:2], *[self.options["md"]["pair_style"],], "2", *pcraw[2:]]) + + lmp.command("pair_style hybrid/scaled v_one %s v_fscale %s"%(self.options["md"]["pair_style"], self.options["md"]["pair_style"])) + lmp.command("pair_coeff %s"%pcnew1) + lmp.command("pair_coeff %s"%pcnew2) + + #start scaling over switching time + lmp.command("fix f1 all nph %s %f %f %f fixedpoint ${xcm} ${ycm} ${zcm}"%(self.iso, pi, + pf, self.options["md"]["pdamp"])) + lmp.command("fix_modify f1 temp tcm") + lmp.command("fix f3 all print 1 \"${dU} $(press) $(vol) ${flambda}\" screen no file forward_%d.dat"%iteration) + lmp.command("run %d"%self.options["md"]["ts"]) + + #unfix + lmp.command("unfix f3") + lmp.command("unfix f1") + + + #switch potential + lmp = ph.set_potential(lmp, self.options) + + #equilibriate scaled hamiltonian + lmp.command("fix f1 all nph %s %f %f %f fixedpoint ${xcm} ${ycm} ${zcm}"%(self.iso, 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") + + + if solid: + solids = ph.find_solid_fraction(os.path.join(self.simfolder, "traj.dat")) + if (solids/lmp.natoms < self.options["conv"]["solid_frac"]): + lmp.close() + raise RuntimeError("System melted, increase size or reduce scaling!") + + + + #reverse scaling + lmp.command("variable flambda equal ramp(${li},${lf})") + lmp.command("variable blambda equal ramp(${lf},${li})") + lmp.command("variable fscale equal v_flambda-1.0") + lmp.command("variable bscale equal v_blambda-1.0") + lmp.command("variable one equal 1.0") + + lmp.command("pair_style hybrid/scaled v_one %s v_bscale %s"%(self.options["md"]["pair_style"], self.options["md"]["pair_style"])) + lmp.command("pair_coeff %s"%pcnew1) + lmp.command("pair_coeff %s"%pcnew2) + + #apply fix and perform switching + lmp.command("fix f1 all nph %s %f %f %f fixedpoint ${xcm} ${ycm} ${zcm}"%(self.iso, pf, + pi, self.options["md"]["pdamp"])) + lmp.command("fix_modify f1 temp tcm") + lmp.command("fix f3 all print 1 \"${dU} $(press) $(vol) ${blambda}\" screen no file backward_%d.dat"%iteration) + lmp.command("run %d"%self.options["md"]["ts"]) + + lmp.command("unfix f3") + lmp.command("unfix f1") + + #close the object + lmp.close() + def integrate_reversible_scaling(self, scale_energy=False, return_values=False): """ Perform integration after reversible scaling diff --git a/calphy/solid.py b/calphy/solid.py index 76eaa71..58468fd 100644 --- a/calphy/solid.py +++ b/calphy/solid.py @@ -446,128 +446,4 @@ def thermodynamic_integration(self): #calculate final free energy self.fe = 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 - 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 up thermostat and barostat and run equilibrium - lmp.command("fix f1 all nph %s %f %f %f"%(self.iso, 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))) - 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 %s %f %f %f fixedpoint ${xcm} ${ycm} ${zcm}"%(self.iso, 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") - - #get output variables - 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 %s %f %f %f fixedpoint ${xcm} ${ycm} ${zcm}"%(self.iso, 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") - - #check for melting - lmp.command("dump 2 all custom 1 traj.dat id type mass x y z vx vy vz") - lmp.command("run 0") - lmp.command("undump 2") - - solids = ph.find_solid_fraction(os.path.join(self.simfolder, "traj.dat")) - if (solids/lmp.natoms < self.options["conv"]["solid_frac"]): - lmp.close() - raise RuntimeError("System melted, increase size or reduce scaling!") - - #equilibriate scaled hamiltonian - lmp.command("fix f1 all nph %s %f %f %f fixedpoint ${xcm} ${ycm} ${zcm}"%( self.iso, 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 %s %f %f %f fixedpoint ${xcm} ${ycm} ${zcm}"%(self.iso, 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() - \ No newline at end of file diff --git a/docs/source/gettingstarted.md b/docs/source/gettingstarted.md index 54bed92..27f96a7 100644 --- a/docs/source/gettingstarted.md +++ b/docs/source/gettingstarted.md @@ -1,44 +1,52 @@ # 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 -In this section, the installation method for `calphy` is explained. Depending on which potential is used, a different version of LAMMPS might be required, leading to a different installation method. The table below lists the pair styles, and the corresponding installation methods. +### Supported operating systems +`calphy` can be installed on Linux and Mac OS based systems. On Windows systems, it is recommended to use [Windows subsystem for Linux](https://docs.microsoft.com/en-us/windows/wsl/install). +### Installation using [Conda](https://anaconda.org/) -| pair style | installation | -|--------------------------------------------------------- |-------------- | -| eam, eam/fs, eam/alloy | [method 1](#id1) | -| pace | [method 2](#id2) | -| adp, snap, sw, tersoff, tersoff/mod, tersoff/modc, meam | [method 3](#id3) | +calphy can be installed directly using [Conda](https://docs.conda.io/en/latest/) from the [conda-forge channel](https://conda-forge.org/) by the following statement: -### Dependencies +``` +conda install -c conda-forge calphy +``` -- lammps 2021.05.27 -- mendeleev 0.7.0 `pip install mendeleev` -- pylammpsmpi 0.0.8 `pip install pylammpsmpi` -- pyscal 2.10.14 `pip install git+https://github.com/srmnitc/pyscal` -- pyyaml 5.4.1 `pip install pyyaml` -- scipy 1.7.0 `pip install scipy` -- tqdm 4.61.2 `pip install tqdm` +### Installation from the repository -#### Optional +calphy can be built from the repository by- -- matplotlib 3.4.2 `pip install matplotlib` -- pytest 6.2.4 `pip install pytest` +``` +git clone https://github.com/ICAMS/calphy.git +cd calphy +python setup.py install --user +``` -### Method 1 +Please note that the list of dependencies mentioned below has to installed manually by the user. -> This method only supports the `eam`, `eam/fs` and `eam/alloy` pair styles. However, direct free energy calculations for solid state can be carried out for any pair style using this installation. +### Using a conda environment It is **strongly** recommended to install and use `calphy` within a conda environment. To see how you can install conda see [here](https://docs.conda.io/projects/conda/en/latest/user-guide/install/). Once a conda distribution is available, the following steps will help set up an environment to use `calphy`. First step is to clone the repository. ``` -git clone https://git.noc.ruhr-uni-bochum.de/atomicclusterexpansion/calphy.git +git clone https://github.com/ICAMS/calphy.git ``` After cloning, an environment can be created from the included file- @@ -48,6 +56,7 @@ cd calphy conda env create -f environment.yml ``` +Note that the conda-forge distribution of LAMMPS will be automatically installed. Alternatively, you can use an existing version of LAMMPS (compiled as library). If a LAMMPS distribution need not be installed, use the `calphy/environment-nolammps.yml` file instead to create the environment. This will install the necessary packages and create an environment called calphy. It can be activated by, ``` @@ -57,44 +66,64 @@ conda activate calphy then, install `calphy` using, ``` -python setup.py install --user +python setup.py install ``` The environment is now set up to run calphy. -### Method 2 +### Dependencies + +- lammps `conda install -c conda-forge lammps` +- mendeleev >=0.7.0 `pip install mendeleev` +- pylammpsmpi >=0.0.8 `pip install pylammpsmpi` +- pyscal >=2.10.14 `pip install git+https://github.com/srmnitc/pyscal` +- pyyaml >=5.4.1 `pip install pyyaml` +- scipy >=1.7.0 `pip install scipy` +- tqdm >=4.61.2 `pip install tqdm` -> This method supports `pace` pair style [1]. Additionally it can also be used if a custom version of LAMMPS is required. Drect free energy calculations for solid state can be carried out for any pair style using this installation. +#### Optional -> **NOTE**: 08 July 2021: The upcoming release of [LAMMPS](https://github.com/lammps/lammps/releases) will significantly change the compilation procedure and render the below information outdated. +- matplotlib >=3.4.2 `pip install matplotlib` +- pytest >=6.2.4 `pip install pytest` +## About [LAMMPS](https://www.lammps.org/) for `calphy` -For using with pair style PACE, the LAMMPS distribution has to be compiled manually: +calphy uses LAMMPS as the driver for molecular dynamics simulations. For calphy to work, LAMMPS needs to be compiled as a library along with the Python interface. The easiest way to do this is to install LAMMPS through the conda-forge channel using: ``` -git clone https://github.com/ICAMS/calphy.git -conda env create -f calphy/environment-nolammps.yml +conda install -c conda-forge lammps ``` -Now proceed with LAMMPS installation +Alternatively, when interatomic potentials with special compilation needs are to be used, LAMMPS (Stable release 29 Sept 2021 and above) can be compiled manually using the following set of instructions: + +First step is to obtain the stable version from [here](https://github.com/lammps/lammps/archive/refs/tags/stable_29Sep2021.tar.gz) and extract the archive. From the extracted archive, the following steps, used in the [conda-forge recipe](https://github.com/conda-forge/lammps-feedstock/blob/master/recipe/build.sh) can be run: ``` -wget https://github.com/lammps/lammps/archive/refs/tags/patch_27May2021.tar.gz -tar -xvf patch_27May2021.tar.gz -cd lammps-patch_27May2021 mkdir build_lib cd build_lib -cmake -D BUILD_LIB=ON -D BUILD_SHARED_LIBS=ON -D BUILD_MPI=ON -D PKG_MANYBODY=ON -D PKG_USER-MISC=ON -D PKG_USER-PACE=ON ../cmake -make -cp liblammps${SHLIB_EXT}* ../src +cmake -D BUILD_LIB=ON -D BUILD_SHARED_LIBS=ON -D BUILD_MPI=ON -D PKG_MPIIO=ON -D LAMMPS_EXCEPTIONS=yes -D PKG_MANYBODY=ON -D PKG_MISC=ON -D PKG_MISC=ON -D PKG_EXTRA-COMPUTE=ON -D PKG_EXTRA-DUMP=ON -D PKG_EXTRA-FIX=ON -D PKG_EXTRA-PAIR=ON ../cmake +make +cp liblammps${SHLIB_EXT}* ../src +cd ../src +make install-python +mkdir -p $PREFIX/include/lammps +cp library.h $PREFIX/include/lammps +cp liblammps${SHLIB_EXT}* "${PREFIX}"/lib/ +cd .. ``` +The above commands only builds the [MANYBODY](https://docs.lammps.org/Packages_details.html#pkg-manybody) package. To use some of the other potentials, the following commands could be added to the `cmake` call. + +- `-D PKG_ML-PACE` for performant [Atomic Cluster Expansion](https://docs.lammps.org/Packages_details.html#pkg-ml-pace) potential. +- `-D PKG_ML-SNAP=ON`for [SNAP potential](https://docs.lammps.org/Packages_details.html#pkg-ml-snap). +- `-D PKG_MEAM=ON` for [MEAM potential](https://docs.lammps.org/Packages_details.html#meam-package). +- `-D PKG_KIM=ON` for [KIM support](https://docs.lammps.org/Packages_details.html#pkg-kim). + Now the python wrapper can be installed: ``` cd ../src -make install-python +make install-python ``` - In the case of a conda environment, the following commands can be used to copy the compiled libraries to an accessible path: ``` @@ -103,7 +132,7 @@ cp library.h $CONDA_PREFIX/include/lammps cp liblammps${SHLIB_EXT}* $CONDA_PREFIX/lib/ ``` -The combined libraries **should be available** on the system or the environment paths. Please see [here](https://lammps.sandia.gov/doc/Python_install.html) for more details. +The combined libraries **should be available** on the system or the environment paths. Please see here for more details. Once LAMMPS is compiled and the libraries are available in an accessible location, the following commands can be used within python to test the installation: @@ -111,96 +140,23 @@ Once LAMMPS is compiled and the libraries are available in an accessible locatio from lammps import lammps lmp = lammps() ``` - -Now, we install `pylammpsmpi` using, +Now, we install pylammpsmpi using, ``` pip install pylammpsmpi ``` - -and finally `calphy`: +and finally calphy: ``` cd calphy python setup.py install --user ``` -### Method 3 - -> This method supports adp, meam, snap, sw, tersoff, tersoff/mod and tersoff/modc pair styles. Additionally it can also be used if a custom version of LAMMPS is required. Drect free energy calculations for solid state can be carried out for any pair style using this installation. - -> **NOTE**: 08 July 2021: The upcoming release of [LAMMPS](https://github.com/lammps/lammps/releases) will significantly change the compilation procedure and render the below information outdated. - -``` -git clone https://github.com/ICAMS/calphy.git -conda env create -f calphy/environment-nolammps.yml -``` - -Now download LAMMPS, - -``` -wget https://github.com/lammps/lammps/archive/refs/tags/patch_27May2021.tar.gz -tar -xvf patch_27May2021.tar.gz -``` - -Get the modified pair styles - -``` -git clone https://github.com/srmnitc/lammps-calphy.git -``` - -The cloned repository contains the necessary files for pair styles. For example, to use sw pair style, replace the original LAMMPS files with the modified ones: - -``` -cp lammps-calphy/src/MANYBODY/pair_sw.* lammps-patch_27May2021/src/MANYBODY/ -``` - -After copying the necessary files, LAMMPS can be installed: - -``` -cd lammps-patch_27May2021 -mkdir build_lib -cd build_lib -cmake -D BUILD_LIB=ON -D BUILD_SHARED_LIBS=ON -D BUILD_MPI=ON -D PKG_MANYBODY=ON -D PKG_USER-MISC=ON ../cmake -make -cp liblammps${SHLIB_EXT}* ../src -``` - -Now the python wrapper can be installed: - -``` -cd ../src -make install-python -``` - -In the case of a conda environment, the following commands can be used to copy the compiled libraries to an accessible path: - -``` -mkdir -p $CONDA_PREFIX/include/lammps -cp library.h $CONDA_PREFIX/include/lammps -cp liblammps${SHLIB_EXT}* $CONDA_PREFIX/lib/ -``` -The combined libraries **should be available** on the system or the environment paths. Please see [here](https://lammps.sandia.gov/doc/Python_install.html) for more details. -Once LAMMPS is compiled and the libraries are available in an accessible location, the following commands can be used within python to test the installation: -``` -from lammps import lammps -lmp = lammps() -``` -Now, we install `pylammpsmpi` using, -``` -pip install pylammpsmpi -``` -and finally `calphy`: -``` -cd calphy -python setup.py install --user -``` -[1] Lysogorskiy, Yury, Cas van der Oord, Anton Bochkarev, Sarath Menon, Matteo Rinaldi, Thomas Hammerschmidt, Matous Mrovec, et al. “Performant Implementation of the Atomic Cluster Expansion (PACE) and Application to Copper and Silicon.” Npj Computational Materials 7, no. 1 (December 2021): 97. https://doi.org/10.1038/s41524-021-00559-9. diff --git a/docs/source/index.rst b/docs/source/index.rst index 0c5723e..fc5f638 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -6,6 +6,9 @@ calphy ====== +calphy is a Python library and command line application which can be used for automated calculation of free energies through molecular dynamics calculations. + + .. toctree:: gettingstarted diff --git a/setup.py b/setup.py index a3f93c8..4e09f77 100644 --- a/setup.py +++ b/setup.py @@ -64,7 +64,7 @@ test_suite='tests', tests_require=test_requirements, url='https://git.noc.ruhr-uni-bochum.de/atomicclusterexpansion/calphy', - version='0.5.8', + version='0.6.0', zip_safe=False, entry_points={ 'console_scripts': [