Calc free energy

`thmd.calc.free_energy.logmfd` ¶

Classes:

LogMFD_FCCUBIC –

Create an Object (class) of Potential, contain some pre-setup information for Energy barrier for LogMFD calculation.

Functions:

replica_logPD_integration –

The function to compute LogPD-based MeanForce
replica_MD_average –

compute Replica_MD_Average from output of MD.
exp_normalize –
read_df –

define lamda(x)
replica_SteerMD –

compute Average Work from output of SteerMD.
find_basin_logmfd_profile –

Args:

`LogMFD_FCCUBIC(Element, Potential_Name, modelType='Bulk', zDirect='001')` ¶

Create an Object (class) of Potential, contain some pre-setup information for Energy barrier for LogMFD calculation. Energy barrier of System does not depend on CV (same value for energy CV: meanCV, nAtomLiquid,...) * Attributes: Element : Cutoff : force cutoff of potential Potential_Name : * Energy barrier Coeff: an estimation of per-atom energy barrier for melting, for a system of N atoms 𝐹=𝑓∗𝑁^(⅔) , then perAtom barrier 𝑓=𝐹/𝑁^(⅔) This method return f(T) * histo_point_coeff: the CV-value at which histogram of solid and liquid meets, coeff of function f(T) = a0 + a1T

Methods: Latt_Const : compute lattice constant at a specific temperature

Element : 'Al', 'Cu',... Potential_Name : 'Cu' : 'Mishin-2001'; 'Foiles-1986';... 'Al' : 'Mishin-1999'; 'Sheng-2011';... modelType : 'Bulk' or 'Plate' zDirect : '001' or '110' or '111' or thickness : thickness of plate

Methods:

meltingBarrier –

for a system of N atoms 𝐹=𝑓∗𝑁^(⅔) , then perAtom barrier 𝑓=𝐹/𝑁^(⅔). This method return f(T)
histo_point –

Value of CV at intersection point of histogram as function y = a + bx

Attributes:

Element –
Potential_Name –
melt_barrier_coeff –
histo_point_coeff –

`Element = Element` `instance-attribute` ¶

`Potential_Name = Potential_Name` `instance-attribute` ¶

`melt_barrier_coeff = melt_barrier_coeff[melt_barrier_key]` `instance-attribute` ¶

`histo_point_coeff = histo_point_coeff` `instance-attribute` ¶

`meltingBarrier(Temp)` ¶

for a system of N atoms 𝐹=𝑓∗𝑁^(⅔) , then perAtom barrier 𝑓=𝐹/𝑁^(⅔). This method return f(T)

`histo_point(Temp)` ¶

Value of CV at intersection point of histogram as function y = a + bx

`replica_logPD_integration(logmfd_files: list[str], replica_files: list[str], beta: float = 1.0)` ¶

The function to compute LogPD-based MeanForce

Parameters:

logmfd_files (list[str]) –

list of logmfd.out files
replica_files (list[str]) –

list of replica.out files
beta (float, default: 1.0 ) –

kB is Boltzmann constant (can be set to 1.0, regardless of kB unit). Defaults to beta = 1.0/(TEMP*kB)

Returns:

file ( obj ) –

contains logPD-based MeanForce

Examples:

free_energy.replica_logPD_intergration(logmfd_files, replica_files)

Requisites:

1. Run logMFD simulations to produce `replica_*/logmfd.out` and `replica_*/replica.out`

```
<logmfd.out>
1:iter_md, 2:Flog(t), …, 6: X(t), 7: V(t), 8: Force(t)
1   F(1), …, X(1), V(1), Force(0)
2   F(3), …, X(2), V(2), Force(1)
```

```
<replica.out>
iter_md, work, weight, cv
1  work(1)   weight(1)  cv(0)
2  work(2)   weight(2)  cv(1)
```

Notes

About the printed values in <replica.out> and <logmfd.out> as emails replied by Tetsuya Morishita. (check thangckt email)

Specify type of function cumulative_trapezoid:np.float64 to be used in numba

cumulative_trapezoid:np.float64 = cumulative_trapezoid,   # this to specify `type` for function `cumulative_trapezoid`
logMFD:pd.DataFrame = read_data.logMFD,
matrix_auto:pd.DataFrame = read_data.matrix_auto):

Refs

[1].https://pubs.acs.org/doi/10.1021/acs.jctc.7b00252 Free Energy Reconstruction from Logarithmic Mean-Force Dynamics Using Multiple Nonequilibrium TrajectoriesFree [2] Exp-normalize trick: https://timvieira.github.io/blog/post/2014/02/11/exp-normalize-trick/

2023 Apr

CV-value in replica.out and logmfd.out are the same, so don't need to shift value in replica.out as below:

Force = [elem['CV1_force'].shift(-1) for elem in logmfd ]

`replica_MD_average(MD_out_files: list[str])` ¶

compute Replica_MD_Average from output of MD.

Parameters:

MD_out_files (list[str]) –

list of "MDout_replica.txt" files

Returns:

–

logPD file: contains logPD-based MeanForce

Requisites

Replica_* files from separate MD simulations

Thang, Jul2020 (update: Sep 2021)

`exp_normalize(x)` ¶

`read_df(file, engine='LAMMPS')` ¶

define lamda(x)

`replica_SteerMD(SteerMD_files, beta=1.0, engine='Lammps')` ¶

compute Average Work from output of SteerMD.

Parameters:

SteerMD_files –

|list| of "SteerMD.txt" files
beta (= 1.0/(TEMP*kB, default: 1.0 ) –

kB is Boltzmann constant (can be set to 1.0, regardless of kB unit)

Returns:

–

aveSteerMD file: contains logPD-based MeanForce

Requisites: 1. Replica_* files from separate MD simulations

Refs

[1]. https://github.com/sandeshkalantre/jarzynski/blob/master/code/Simulations%20on%20Harmonic%20Oscillator%20Model.ipynb [2]. https://www.plumed.org/doc-v2.6/user-doc/html/belfast-5.html#belfast-5-work [3]. Exp-normalize trick: https://timvieira.github.io/blog/post/2014/02/11/exp-normalize-trick/

`find_basin_logmfd_profile(ldf: list[pd.DataFrame], mfd_interval: list[float], hill_region: list = [0.2, 0.6]) -> list[pd.DataFrame]` ¶

Parameters:

ldf (list[DataFrame]) –

list of pd.DataFrames of LogMFD.out files
mfd_interval (list) –

interval of MFDstep in PLUMED input files
hill_region (1x2 list, default: [0.2, 0.6] ) –

the region of CV that contains saddle point. The Left/Right basin is located at left/right of this region. This region should exclude the bassins.

Returns:

df_min_right ( DataFrame ) –

DataFrame of right basin minima
df_min_left ( DataFrame ) –

DataFrame of left basin minima
df_max ( DataFrame ) –

DataFrame of maxima energy
df_diff ( DataFrame ) –

DataFrame of energy difference between minima and maxima

`thmd.calc.free_energy.TImethod` ¶

Functions:

AS_integration –

Compute Free energy difference in TI method - Lambda integration along Adiabatic Switching path
RS_integration –

Compute Free energy difference in TI method - Lambda integration along Reversible Scaling path
Helmholtz_excess_UF –

this func. is to compute the excess Helmholtz freeEnergy as Eq.(25) by R.Paula Leite 2016.

Attributes:

kB –

`kB = sc.value('Boltzmann constant in eV/K')` `module-attribute` ¶

`AS_integration(forwFile, backFile)` ¶

Compute Free energy difference in TI method - Lambda integration along Adiabatic Switching path

Parameters:

forwFile (str) –

file name of forward integration (must contain 2 first columns are: dE, lambda)
backFile (str) –

file name of backward integration (must contain 2 first columns are: dE, lambda)

Returns:

W ( float ) –

irreversible TI work
Q ( float ) –

dissipation

`RS_integration(forwFile, backFile, T0=0, F0=0, kB=kB)` ¶

Compute Free energy difference in TI method - Lambda integration along Reversible Scaling path

Parameters:

forwFile (str) –

file name of forward integration (must contain 2 first columns are: dE, lambda)
backFile (str) –

file name of backward integration (must contain 2 first columns are: dE, lambda)
T0 (float, default: 0 ) –

reference temperature. Default is 0.
F0 (float, default: 0 ) –

free energy at reference temperature. Default is 0.
kB (float, default: kB ) –

Boltzmann constant (chose unit to consist with F0). Default is kB in unit eV/K.

Returns:

T ( array ) –

temperature
Ft ( array ) –

free energy as a function of temperature.
W ( array ) –

cumulative work

`Helmholtz_excess_UF(p, x)` ¶

this func. is to compute the excess Helmholtz freeEnergy as Eq.(25) by R.Paula Leite 2016. J.Chem.Phys.145,no.19:194101. https://doi.org/10.1063/1.4967775. p : # UFM p-parameter x = b*rho : # adimensional variable, involved UFM-sigma parameter, rho in unit [1/A^3]

the excess_Helmholtz free enegy in eV : (beta*Fexc)/N

the pressure : beta*b*P

Calc free energy

thmd.calc.free_energy.logmfd ¶

LogMFD_FCCUBIC(Element, Potential_Name, modelType='Bulk', zDirect='001') ¶

Element = Element instance-attribute ¶

Potential_Name = Potential_Name instance-attribute ¶

melt_barrier_coeff = melt_barrier_coeff[melt_barrier_key] instance-attribute ¶

histo_point_coeff = histo_point_coeff instance-attribute ¶

meltingBarrier(Temp) ¶

histo_point(Temp) ¶

replica_logPD_integration(logmfd_files: list[str], replica_files: list[str], beta: float = 1.0) ¶

replica_MD_average(MD_out_files: list[str]) ¶

exp_normalize(x) ¶

read_df(file, engine='LAMMPS') ¶

replica_SteerMD(SteerMD_files, beta=1.0, engine='Lammps') ¶

find_basin_logmfd_profile(ldf: list[pd.DataFrame], mfd_interval: list[float], hill_region: list = [0.2, 0.6]) -> list[pd.DataFrame] ¶

thmd.calc.free_energy.TImethod ¶

kB = sc.value('Boltzmann constant in eV/K') module-attribute ¶

AS_integration(forwFile, backFile) ¶

RS_integration(forwFile, backFile, T0=0, F0=0, kB=kB) ¶

Helmholtz_excess_UF(p, x) ¶

`thmd.calc.free_energy.logmfd` ¶

`LogMFD_FCCUBIC(Element, Potential_Name, modelType='Bulk', zDirect='001')` ¶

`Element = Element` `instance-attribute` ¶

`Potential_Name = Potential_Name` `instance-attribute` ¶

`melt_barrier_coeff = melt_barrier_coeff[melt_barrier_key]` `instance-attribute` ¶

`histo_point_coeff = histo_point_coeff` `instance-attribute` ¶

`meltingBarrier(Temp)` ¶

`histo_point(Temp)` ¶

`replica_logPD_integration(logmfd_files: list[str], replica_files: list[str], beta: float = 1.0)` ¶

`replica_MD_average(MD_out_files: list[str])` ¶

`exp_normalize(x)` ¶

`read_df(file, engine='LAMMPS')` ¶

`replica_SteerMD(SteerMD_files, beta=1.0, engine='Lammps')` ¶

`find_basin_logmfd_profile(ldf: list[pd.DataFrame], mfd_interval: list[float], hill_region: list = [0.2, 0.6]) -> list[pd.DataFrame]` ¶

`thmd.calc.free_energy.TImethod` ¶

`kB = sc.value('Boltzmann constant in eV/K')` `module-attribute` ¶

`AS_integration(forwFile, backFile)` ¶

`RS_integration(forwFile, backFile, T0=0, F0=0, kB=kB)` ¶

`Helmholtz_excess_UF(p, x)` ¶