2020import ase .io
2121import chemiscope
2222import ipi
23- import ipi .utils .parsing as pimdmooc
2423import matplotlib .pyplot as plt
2524import numpy as np
2625
@@ -110,7 +109,8 @@ def PES(x, y, z):
110109# The harmonic approximation to the PES is essentially
111110# a truncated Taylor series expansion to second
112111# order around one of its minima.
113- # :math:`V^{\text{harm}} = V(q_0) + \frac{2}{2} \left.\frac{\partial^2 V}{\partial q^2}\right|_{q_0} (q - q_0)^2`
112+ # :math:`V^{\text{harm}} = V(q_0) +
113+ # \frac{1}{2} \left.\frac{\partial^2 V}{\partial q^2}\right|_{q_0} (q - q_0)^2`
114114# where :math:`q = (x,y,z)`
115115# is a position vector and :math:`q_0 = \text{arg min}_{q} V` is
116116# the position where the PES has a local minimum.
@@ -150,7 +150,7 @@ def PES(x, y, z):
150150# simply move towards a local minimum of the PES. There are
151151# many algorithms for locally optimizing high-dimensional functions;
152152# here we use the robust
153- # `BFGS <https://en.wikipedia.org/wiki/ Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm>`_
153+ # `BFGS` ( Broyden-Fletcher-Goldfarb-Shanno) quasi-Newton
154154# method. The tolerances set thershold values for changes in the energy,
155155# positions and forces, that are sufficient to deem an optimization converged.
156156#
@@ -252,8 +252,10 @@ def PES(x, y, z):
252252# - \epsilon} \right)
253253#
254254# At every step instead of performing "dynamics", we will displace a degree of
255- # freedom along :math:`\pm \epsilon` and estimate one row of the Hessian. `pos_shift` sets
256- # the value of :math:`\epsilon` while `asr` zeros out the blocks of the Hessian due to continuous
255+ # freedom along :math:`\pm \epsilon` and estimate one row of the Hessian.
256+ # `pos_shift` sets
257+ # the value of :math:`\epsilon` while `asr` zeros out the blocks
258+ # of the Hessian due to continuous
257259# symmetries (translations or rotations for solids or clusters). In this example,
258260# we set this option to `none` as our system doesn't possess any continuous symmetries.
259261#
@@ -396,9 +398,11 @@ def quantum_harmonic_free_energy(Ws, T):
396398#
397399# Calculating free energies beyond the harmonic approximation is non-trivial.
398400# There exist a familty of methods that can solve the vibrational Schroedinger
399- # Equation by approximating the anharmonic component of the PES, yielding an amharmonic
401+ # Equation by approximating the anharmonic component of
402+ # the PES, yielding an amharmonic
400403# free energy. While highly effective for low-dimensional or mildly anharmonic systems,
401- # the method of resort for *numerically-exact amharmonic free energies* of solid and clusters
404+ # the method of resort for *numerically-exact amharmonicfree energies*
405+ # of solid and clusters
402406# is the thermodynamic integration method combined with the path-integral method
403407# ( for applications see Refs.
404408# `M. Rossi et al, PRL (2016) <https://doi.org/10.1103/PhysRevLett.117.115702>`_,
@@ -439,7 +443,8 @@ def quantum_harmonic_free_energy(Ws, T):
439443# Classical Statistical Mechanics
440444# -------------------------------
441445#
442- # Let's first compute the anharmonic free energy difference within the classical approximation.
446+ # Let's first compute the anharmonic free
447+ # energy difference within the classical approximation.
443448#
444449# To create the input file for I-PI we need to defines a "mixed"
445450# :math:`\lambda`-dependent Hamiltonian. Let's see the new most important parts.
@@ -463,8 +468,12 @@ def quantum_harmonic_free_energy(Ws, T):
463468# <!--> defines the harmonic PES <-->
464469# <ffdebye name='debye'>
465470# <hessian shape='(3,3)' mode='file'> HESSIAN_FILE </hessian>
466- # <x_reference mode='file'> X_REF_FILE <!--> relative path to a file containing the optimized positon vector <--> </x_reference>
467- # <v_reference> V_REF <!--> the value of the PES at its local minimum <--> </v_reference>
471+ # <x_reference mode='file'> X_REF_FILE
472+ # <!--> relative path to a file containing the optimized positon vector <-->
473+ # </x_reference>
474+ # <v_reference> V_REF
475+ # <!--> the value of the PES at its local minimum <-->
476+ # </v_reference>
468477# </ffdebye>
469478#
470479# an intrinsic harmonic forcefield that builds the harmonic potential.
@@ -494,10 +503,10 @@ def quantum_harmonic_free_energy(Ws, T):
494503#
495504# .. code-block:: xml
496505#
497- # <forces>
498- # <force forcefield='debye' weight=''> </force> <!--> set this to lambda <-->
499- # <force forcefield='driver' weight=''> </force> <!--> set this to 1 - lambda <-->
500- # </forces>
506+ # <forces>
507+ # <force forcefield='debye' weight=''> </force> <!--> set this to lambda <-->
508+ # <force forcefield='driver' weight=''> </force> <!--> set this to 1 - lambda <-->
509+ # </forces>
501510#
502511# You can print out the harmonic and the anharmonic
503512# components as a <property> in the output class
@@ -521,9 +530,14 @@ def quantum_harmonic_free_energy(Ws, T):
521530 with open (f"cti/input_{ i } .xml" , "w" ) as f :
522531 strfile = """<simulation verbosity='low'>
523532 <output prefix='cti_{i}'>
524- <properties filename='out' stride='10' flush='10'> [ step, time{{picosecond}}, conserved, temperature{{kelvin}}, kinetic_cv, potential, pressure_cv, volume, ensemble_temperature ] </properties>
525- <properties filename='pots' stride='10' flush='10'> [ pot_component_raw(0), pot_component_raw(1) ] </properties>
526- <trajectory filename='pos1' stride='10' bead='0' flush='10'> positions </trajectory>
533+ <properties filename='out' stride='10' flush='10'> [ step, time{{picosecond}},
534+ conserved, temperature{{kelvin}}, kinetic_cv, potential, pressure_cv,
535+ volume, ensemble_temperature ] </properties>
536+ <properties filename='pots' stride='10' flush='10'>
537+ [ pot_component_raw(0), pot_component_raw(1) ] </properties>
538+ <trajectory filename='pos1' stride='10' bead='0' flush='10'>
539+ positions
540+ </trajectory>
527541 <trajectory filename='xc' stride='10' flush='10'> x_centroid </trajectory>
528542 <checkpoint stride='4000'/>
529543 </output>
@@ -534,21 +548,21 @@ def quantum_harmonic_free_energy(Ws, T):
534548 <latency> 1e-3 </latency>
535549 </ffsocket>
536550 <ffdebye name='debye'>
537- <hessian shape='(3,3)' mode='file'> {HESSIAN_FILE} </hessian>
538- <x_reference mode='file'> {X_REF_FILE} </x_reference>
539- <v_reference> {V_REF} </v_reference>
551+ <hessian shape='(3,3)' mode='file'> {HESSIAN_FILE} </hessian>
552+ <x_reference mode='file'> {X_REF_FILE} </x_reference>
553+ <v_reference> {V_REF} </v_reference>
540554 </ffdebye>
541555 <system>
542556 <initialize nbeads='1'>
543- <file mode='xyz'> ./cti/init.xyz </file>
557+ <file mode='xyz'> ./cti/init.xyz </file>
544558 <velocities mode='thermal' units='kelvin'> 300 </velocities>
545559 </initialize>
546560 <forces>
547561 <force forcefield='debye' weight='{i}'> </force>
548562 <force forcefield='driver' weight='{im1:2.1f}'> </force>
549563 </forces>
550564 <motion mode='dynamics'>
551- <fixcom> False </fixcom>
565+ <fixcom> False </fixcom>
552566 <dynamics mode='nvt'>
553567 <timestep units='femtosecond'> 1.00 </timestep>
554568 <thermostat mode='pile_l'>
@@ -584,7 +598,7 @@ def quantum_harmonic_free_energy(Ws, T):
584598#
585599# wait 5
586600#
587- # for x in {0..10..2 }; do
601+ # for x in {0..10}; do
588602# i-pi-driver -u -h f${x} -m doublewell_1D &
589603# done
590604#
@@ -613,7 +627,7 @@ def quantum_harmonic_free_energy(Ws, T):
613627duerr_list = []
614628
615629dir_list = ["0.0" , "0.2" , "0.4" , "0.6" , "0.8" , "1.0" ]
616- l_list = [float (l ) for l in dir_list ]
630+ l_list = [float (lst ) for lst in dir_list ]
617631
618632for i in dir_list :
619633
@@ -702,9 +716,14 @@ def classical_harmonic_free_energy(Ws, T):
702716 with open (f"qti/input_{ i } .xml" , "w" ) as f :
703717 strfile = """<simulation verbosity='low'>
704718 <output prefix='qti_{i}'>
705- <properties filename='out' stride='10' flush='10'> [ step, time{{picosecond}}, conserved, temperature{{kelvin}}, kinetic_cv, potential, pressure_cv, volume, ensemble_temperature ] </properties>
706- <properties filename='pots' stride='10' flush='10'> [ pot_component_raw(0), pot_component_raw(1) ] </properties>
707- <trajectory filename='pos1' stride='10' bead='0' flush='10'> positions </trajectory>
719+ <properties filename='out' stride='10' flush='10'> [ step, time{{picosecond}},
720+ conserved, temperature{{kelvin}}, kinetic_cv, potential,
721+ pressure_cv, volume, ensemble_temperature ] </properties>
722+ <properties filename='pots' stride='10' flush='10'>
723+ [ pot_component_raw(0), pot_component_raw(1) ] </properties>
724+ <trajectory filename='pos1' stride='10' bead='0' flush='10'>
725+ positions
726+ </trajectory>
708727 <trajectory filename='xc' stride='10' flush='10'> x_centroid </trajectory>
709728 <checkpoint stride='4000'/>
710729 </output>
@@ -715,21 +734,21 @@ def classical_harmonic_free_energy(Ws, T):
715734 <latency> 1e-3 </latency>
716735 </ffsocket>
717736 <ffdebye name='debye'>
718- <hessian shape='(3,3)' mode='file'> {HESSIAN_FILE} </hessian>
719- <x_reference mode='file'> {X_REF_FILE} </x_reference>
720- <v_reference> {V_REF} </v_reference>
737+ <hessian shape='(3,3)' mode='file'> {HESSIAN_FILE} </hessian>
738+ <x_reference mode='file'> {X_REF_FILE} </x_reference>
739+ <v_reference> {V_REF} </v_reference>
721740 </ffdebye>
722741 <system>
723742 <initialize nbeads='32'>
724- <file mode='xyz'> ./qti/init.xyz </file>
743+ <file mode='xyz'> ./qti/init.xyz </file>
725744 <velocities mode='thermal' units='kelvin'> 300 </velocities>
726745 </initialize>
727746 <forces>
728747 <force forcefield='debye' weight='{i}'> </force>
729748 <force forcefield='driver' weight='{im1:2.1f}'> </force>
730749 </forces>
731750 <motion mode='dynamics'>
732- <fixcom> False </fixcom>
751+ <fixcom> False </fixcom>
733752 <dynamics mode='nvt'>
734753 <timestep units='femtosecond'> 1.00 </timestep>
735754 <thermostat mode='pile_l'>
@@ -781,7 +800,7 @@ def classical_harmonic_free_energy(Ws, T):
781800Q_duerr_list = []
782801
783802Q_dir_list = ["0.0" , "0.2" , "0.4" , "0.6" , "0.8" , "1.0" ]
784- Q_l_list = [float (l ) for l in Q_dir_list ]
803+ Q_l_list = [float (lst ) for lst in Q_dir_list ]
785804
786805for i in Q_dir_list :
787806 filename = f"qti_{ i } .pots"
0 commit comments