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

    '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

1. About the printed values in <replica.out> and <logmfd.out> as emails replied by Tetsuya Morishita. (check thangckt email)
2. 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

1. 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) –

  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

Used by:

• thmd calc Calc free energy  TImethod  RS_integration

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]

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

1. the pressure : beta*b*P