20
20
import ase .io
21
21
import chemiscope
22
22
import ipi
23
- import ipi .utils .parsing as pimdmooc
24
23
import matplotlib .pyplot as plt
25
24
import numpy as np
26
25
@@ -110,7 +109,8 @@ def PES(x, y, z):
110
109
# The harmonic approximation to the PES is essentially
111
110
# a truncated Taylor series expansion to second
112
111
# 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`
114
114
# where :math:`q = (x,y,z)`
115
115
# is a position vector and :math:`q_0 = \text{arg min}_{q} V` is
116
116
# the position where the PES has a local minimum.
@@ -150,7 +150,7 @@ def PES(x, y, z):
150
150
# simply move towards a local minimum of the PES. There are
151
151
# many algorithms for locally optimizing high-dimensional functions;
152
152
# 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
154
154
# method. The tolerances set thershold values for changes in the energy,
155
155
# positions and forces, that are sufficient to deem an optimization converged.
156
156
#
@@ -252,8 +252,10 @@ def PES(x, y, z):
252
252
# - \epsilon} \right)
253
253
#
254
254
# 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
257
259
# symmetries (translations or rotations for solids or clusters). In this example,
258
260
# we set this option to `none` as our system doesn't possess any continuous symmetries.
259
261
#
@@ -396,9 +398,11 @@ def quantum_harmonic_free_energy(Ws, T):
396
398
#
397
399
# Calculating free energies beyond the harmonic approximation is non-trivial.
398
400
# 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
400
403
# 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
402
406
# is the thermodynamic integration method combined with the path-integral method
403
407
# ( for applications see Refs.
404
408
# `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):
439
443
# Classical Statistical Mechanics
440
444
# -------------------------------
441
445
#
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.
443
448
#
444
449
# To create the input file for I-PI we need to defines a "mixed"
445
450
# :math:`\lambda`-dependent Hamiltonian. Let's see the new most important parts.
@@ -463,8 +468,12 @@ def quantum_harmonic_free_energy(Ws, T):
463
468
# <!--> defines the harmonic PES <-->
464
469
# <ffdebye name='debye'>
465
470
# <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>
468
477
# </ffdebye>
469
478
#
470
479
# an intrinsic harmonic forcefield that builds the harmonic potential.
@@ -494,10 +503,10 @@ def quantum_harmonic_free_energy(Ws, T):
494
503
#
495
504
# .. code-block:: xml
496
505
#
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>
501
510
#
502
511
# You can print out the harmonic and the anharmonic
503
512
# components as a <property> in the output class
@@ -521,9 +530,14 @@ def quantum_harmonic_free_energy(Ws, T):
521
530
with open (f"cti/input_{ i } .xml" , "w" ) as f :
522
531
strfile = """<simulation verbosity='low'>
523
532
<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>
527
541
<trajectory filename='xc' stride='10' flush='10'> x_centroid </trajectory>
528
542
<checkpoint stride='4000'/>
529
543
</output>
@@ -534,21 +548,21 @@ def quantum_harmonic_free_energy(Ws, T):
534
548
<latency> 1e-3 </latency>
535
549
</ffsocket>
536
550
<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>
540
554
</ffdebye>
541
555
<system>
542
556
<initialize nbeads='1'>
543
- <file mode='xyz'> ./cti/init.xyz </file>
557
+ <file mode='xyz'> ./cti/init.xyz </file>
544
558
<velocities mode='thermal' units='kelvin'> 300 </velocities>
545
559
</initialize>
546
560
<forces>
547
561
<force forcefield='debye' weight='{i}'> </force>
548
562
<force forcefield='driver' weight='{im1:2.1f}'> </force>
549
563
</forces>
550
564
<motion mode='dynamics'>
551
- <fixcom> False </fixcom>
565
+ <fixcom> False </fixcom>
552
566
<dynamics mode='nvt'>
553
567
<timestep units='femtosecond'> 1.00 </timestep>
554
568
<thermostat mode='pile_l'>
@@ -584,7 +598,7 @@ def quantum_harmonic_free_energy(Ws, T):
584
598
#
585
599
# wait 5
586
600
#
587
- # for x in {0..10..2 }; do
601
+ # for x in {0..10}; do
588
602
# i-pi-driver -u -h f${x} -m doublewell_1D &
589
603
# done
590
604
#
@@ -613,7 +627,7 @@ def quantum_harmonic_free_energy(Ws, T):
613
627
duerr_list = []
614
628
615
629
dir_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 ]
617
631
618
632
for i in dir_list :
619
633
@@ -702,9 +716,14 @@ def classical_harmonic_free_energy(Ws, T):
702
716
with open (f"qti/input_{ i } .xml" , "w" ) as f :
703
717
strfile = """<simulation verbosity='low'>
704
718
<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>
708
727
<trajectory filename='xc' stride='10' flush='10'> x_centroid </trajectory>
709
728
<checkpoint stride='4000'/>
710
729
</output>
@@ -715,21 +734,21 @@ def classical_harmonic_free_energy(Ws, T):
715
734
<latency> 1e-3 </latency>
716
735
</ffsocket>
717
736
<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>
721
740
</ffdebye>
722
741
<system>
723
742
<initialize nbeads='32'>
724
- <file mode='xyz'> ./qti/init.xyz </file>
743
+ <file mode='xyz'> ./qti/init.xyz </file>
725
744
<velocities mode='thermal' units='kelvin'> 300 </velocities>
726
745
</initialize>
727
746
<forces>
728
747
<force forcefield='debye' weight='{i}'> </force>
729
748
<force forcefield='driver' weight='{im1:2.1f}'> </force>
730
749
</forces>
731
750
<motion mode='dynamics'>
732
- <fixcom> False </fixcom>
751
+ <fixcom> False </fixcom>
733
752
<dynamics mode='nvt'>
734
753
<timestep units='femtosecond'> 1.00 </timestep>
735
754
<thermostat mode='pile_l'>
@@ -781,7 +800,7 @@ def classical_harmonic_free_energy(Ws, T):
781
800
Q_duerr_list = []
782
801
783
802
Q_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 ]
785
804
786
805
for i in Q_dir_list :
787
806
filename = f"qti_{ i } .pots"
0 commit comments