rename example folders

This commit is contained in:
Axel Kohlmeyer
2021-06-29 11:23:47 -04:00
parent b8c5f08e75
commit 4342987043
1351 changed files with 4 additions and 3 deletions

View File

@ -0,0 +1,10 @@
There are example scripts for using the USER-QTB package.
See its src/USER-QTB/README file for more info on
quantum nuclear effects and when they are important
to include in your model.
Authors Information and Contact
The person who created this package and these examples is Yuan Shen
(sy0302 at stanford.edu) at Stanford University. Contact him directly
if you have questions.

View File

@ -0,0 +1,15 @@
#This script implements the BKS pair potential for various silicon dioxide compounds. Inner part is fixed with a harmonic potential. Long range Coulomb interactions are evaluated with the pppm method.
#Pair Potentials
pair_style hybrid/overlay buck/coul/long ${cut_off} table linear 39901
pair_coeff 1 1 buck/coul/long 0.0 1.0 0.0 #No interactions between Si atoms
pair_coeff 1 2 buck/coul/long 18003.757200 0.205205 133.538100
pair_coeff 2 2 buck/coul/long 1388.773000 0.362319 175.000000 #BKS interaction in PRL 64 1955 (1990)
pair_modify shift yes
pair_coeff 1 2 table potential_SiO2.TPF Si-O ${cut_off}
pair_coeff 2 2 table potential_SiO2.TPF O-O ${cut_off} #See the potential file for more information
kspace_style pppm 1.0e-4
#Neighbor style
neighbor 2.0 bin
neigh_modify check yes every 1 delay 0 page 100000 one 2000

View File

@ -0,0 +1,62 @@
## This script first constructs an alpha quartz structure of a given size. It then uses fix qtb to equilibrate the computational cell to the specified temperature and pressure.
## This part defines units, alpha-quartz crystal, and atomic information
#General
units metal
dimension 3
boundary p p p
atom_style charge
#Lattice
lattice custom 1.0 &
a1 4.916000 0.000000 0.000000 &
a2 -2.45800 4.257381 0.000000 &
a3 0.000000 0.000000 5.405400 &
&
basis 0.469700 0.000000 0.000000 &
basis 0.000000 0.469700 0.666667 &
basis 0.530300 0.530300 0.333333 &
&
basis 0.413500 0.266900 0.119100 &
basis 0.266900 0.413500 0.547567 &
basis 0.733100 0.146600 0.785767 &
basis 0.586500 0.853400 0.214233 &
basis 0.853400 0.586500 0.452433 &
basis 0.146600 0.733100 0.880900 #American Mineralogist 65 920 1980 (Space Group 154)
#Computational Cell
region orthorhombic_unit_cell block 0 4.916000 0 8.514762 0 5.405400 units box
create_box 2 orthorhombic_unit_cell
create_atoms 1 box &
basis 1 1 &
basis 2 1 &
basis 3 1 &
basis 4 2 &
basis 5 2 &
basis 6 2 &
basis 7 2 &
basis 8 2 &
basis 9 2
replicate ${x_rep} ${y_rep} ${z_rep}
#Atomic Information
mass 1 28.085500
mass 2 15.999400
set type 1 charge +2.4
set type 2 charge -1.2
## This part implements the BKS pair potential with a cut-off distance for the Buckingham term. Long range Coulomb interactions are evaluated with the pppm method.
include alpha_quartz_potential.mod
## This part equilibrates your crystal to a pressure of ${pressure}(unit pressure) and a temperature of ${temperature}(unit temperatureture) with quantum nuclear effects
variable p_damp equal ${delta_t}*1000 #Recommended pressure damping parameter in fix nph
fix scapegoat_qtb all nph iso ${pressure} ${pressure} ${p_damp} ptemp ${temperature} #NPH does the time integration
fix quartz_qtb all qtb temp ${temperature} damp ${damp_qtb} seed 35082 f_max 120.00 N_f 100 #Change f_max (THz) if your Debye frequency is higher
thermo_style custom step temp press etotal vol lx ly lz pxx pyy pzz pxy pyz pxz
thermo 200
run 2000 # 2 ps
unfix quartz_qtb
unfix scapegoat_qtb

View File

@ -0,0 +1,54 @@
## This script first uses fix qtb to equilibrate alpha quartz structure to an initial state with quantum nuclear correction and then simulate shock induced phase transition through the quantum thermal bath multi-scale shock technique
variable x_rep equal 2 #plot is made with x_rep = 8 #x-direction replication number
variable y_rep equal 1 #plot is made with y_rep = 5 #y-direction replication number
variable z_rep equal 4 #plot is made with z_rep = 15 #z-direction replication number
variable cut_off equal 10.0 #Cut-off distance for the Buckingham term (Angstrom in metal units)
variable pressure equal 1.03125 #Initial state pressure (bar in metal units)
variable temperature equal 300.0 #Initial state quantum temperature (K in metal units)
variable delta_t equal 1.0e-3 #MD timestep length (ps in metal units)
variable damp_qtb equal 1 #1/gamma where gamma is the friction coefficient in quantum thermal bath (ps in metal units)
variable v_msst equal 78.0 #Shock velocity (Angstrom/ps in metal units)
variable q_msst equal 40.0 #Box mass-like parameter in the MSST (mass^2/length^4, where mass=grams/mole and length=Angstrom in metal units)
variable tscale_msst equal 0.05 #Temperature reduction parameter in the MSST (unitless)
variable eta_qbmsst equal 1.0 #Coupling constant between the shock and the quantum thermal bath (unitless constant)
## This part uses fix qtb to prepare alpha-quartz with quantum nuclear correction of the initial state
include alpha_quartz_qtb.mod
## This part demonstrates how to retart fix qbmsst during any stage of the shock simulation.
## PPPM may break down when compression ratio becomes extremely large. One can always use this restart technique to resume the shock simulation.
#Compression restart 1
reset_timestep 0
#Beta is the number of time steps between each update of the quantum bath temperature. Setting a larger beta can reduce thermal flactuations.
fix shock all qbmsst z ${v_msst} q ${q_msst} tscale ${tscale_msst} damp ${damp_qtb} f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
variable dhug equal f_shock[1]
variable dray equal f_shock[2]
variable lgr_vel equal f_shock[3]
variable lgr_pos equal f_shock[4]
variable T_qm equal f_shock[5] #Temperature with quantum nuclear correction
thermo_style custom step v_T_qm press econserve vol lx ly lz pzz v_dhug v_dray v_lgr_vel v_lgr_pos
thermo 200
timestep ${delta_t}
run 1000
write_restart restart.1000
#Compression restart 2
#Read restart file and load potential again
clear
read_restart restart.1000
include alpha_quartz_potential.mod
#Use the same fix id and add no tscale if the system is already compressed
fix shock all qbmsst z ${v_msst} q ${q_msst} tscale 0.0 damp ${damp_qtb} f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
variable dhug equal f_shock[1]
variable dray equal f_shock[2]
variable lgr_vel equal f_shock[3]
variable lgr_pos equal f_shock[4]
variable T_qm equal f_shock[5] #Temperature with quantum nuclear correction
thermo_style custom step v_T_qm press econserve vol lx ly lz pzz v_dhug v_dray v_lgr_vel v_lgr_pos
thermo 500
timestep ${delta_t}
#restart 1000 restart
run 10000 #10 ps
shell rm restart.1000

View File

@ -0,0 +1,435 @@
LAMMPS (24 Dec 2020)
using 1 OpenMP thread(s) per MPI task
## This script first uses fix qtb to equilibrate alpha quartz structure to an initial state with quantum nuclear correction and then simulate shock induced phase transition through the quantum thermal bath multi-scale shock technique
variable x_rep equal 2 #plot is made with x_rep = 8 #x-direction replication number
variable y_rep equal 1 #plot is made with y_rep = 5 #y-direction replication number
variable z_rep equal 4 #plot is made with z_rep = 15 #z-direction replication number
variable cut_off equal 10.0 #Cut-off distance for the Buckingham term (Angstrom in metal units)
variable pressure equal 1.03125 #Initial state pressure (bar in metal units)
variable temperature equal 300.0 #Initial state quantum temperature (K in metal units)
variable delta_t equal 1.0e-3 #MD timestep length (ps in metal units)
variable damp_qtb equal 1 #1/gamma where gamma is the friction coefficient in quantum thermal bath (ps in metal units)
variable v_msst equal 78.0 #Shock velocity (Angstrom/ps in metal units)
variable q_msst equal 40.0 #Box mass-like parameter in the MSST (mass^2/length^4, where mass=grams/mole and length=Angstrom in metal units)
variable tscale_msst equal 0.05 #Temperature reduction parameter in the MSST (unitless)
variable eta_qbmsst equal 1.0 #Coupling constant between the shock and the quantum thermal bath (unitless constant)
## This part uses fix qtb to prepare alpha-quartz with quantum nuclear correction of the initial state
include alpha_quartz_qtb.mod
## This script first constructs an alpha quartz structure of a given size. It then uses fix qtb to equilibrate the computational cell to the specified temperature and pressure.
## This part defines units, alpha-quartz crystal, and atomic information
#General
units metal
dimension 3
boundary p p p
atom_style charge
#Lattice
lattice custom 1.0 a1 4.916000 0.000000 0.000000 a2 -2.45800 4.257381 0.000000 a3 0.000000 0.000000 5.405400 basis 0.469700 0.000000 0.000000 basis 0.000000 0.469700 0.666667 basis 0.530300 0.530300 0.333333 basis 0.413500 0.266900 0.119100 basis 0.266900 0.413500 0.547567 basis 0.733100 0.146600 0.785767 basis 0.586500 0.853400 0.214233 basis 0.853400 0.586500 0.452433 basis 0.146600 0.733100 0.880900 #American Mineralogist 65 920 1980 (Space Group 154)
Lattice spacing in x,y,z = 7.3740000 4.2573810 5.4054000
#Computational Cell
region orthorhombic_unit_cell block 0 4.916000 0 8.514762 0 5.405400 units box
create_box 2 orthorhombic_unit_cell
Created orthogonal box = (0.0000000 0.0000000 0.0000000) to (4.9160000 8.5147620 5.4054000)
1 by 1 by 1 MPI processor grid
create_atoms 1 box basis 1 1 basis 2 1 basis 3 1 basis 4 2 basis 5 2 basis 6 2 basis 7 2 basis 8 2 basis 9 2
Created 18 atoms
create_atoms CPU = 0.000 seconds
replicate ${x_rep} ${y_rep} ${z_rep}
replicate 2 ${y_rep} ${z_rep}
replicate 2 1 ${z_rep}
replicate 2 1 4
Replicating atoms ...
orthogonal box = (0.0000000 0.0000000 0.0000000) to (9.8320000 8.5147620 21.621600)
1 by 1 by 1 MPI processor grid
144 atoms
replicate CPU = 0.001 seconds
#Atomic Information
mass 1 28.085500
mass 2 15.999400
set type 1 charge +2.4
Setting atom values ...
48 settings made for charge
set type 2 charge -1.2
Setting atom values ...
96 settings made for charge
## This part implements the BKS pair potential with a cut-off distance for the Buckingham term. Long range Coulomb interactions are evaluated with the pppm method.
include alpha_quartz_potential.mod
#This script implements the BKS pair potential for various silicon dioxide compounds. Inner part is fixed with a harmonic potential. Long range Coulomb interactions are evaluated with the pppm method.
#Pair Potentials
pair_style hybrid/overlay buck/coul/long ${cut_off} table linear 39901
pair_style hybrid/overlay buck/coul/long 10 table linear 39901
pair_coeff 1 1 buck/coul/long 0.0 1.0 0.0 #No interactions between Si atoms
pair_coeff 1 2 buck/coul/long 18003.757200 0.205205 133.538100
pair_coeff 2 2 buck/coul/long 1388.773000 0.362319 175.000000 #BKS interaction in PRL 64 1955 (1990)
pair_modify shift yes
pair_coeff 1 2 table potential_SiO2.TPF Si-O ${cut_off}
pair_coeff 1 2 table potential_SiO2.TPF Si-O 10
pair_coeff 2 2 table potential_SiO2.TPF O-O ${cut_off} #See the potential file for more information
pair_coeff 2 2 table potential_SiO2.TPF O-O 10
WARNING: 1 of 39901 force values in table O-O are inconsistent with -dE/dr.
Should only be flagged at inflection points (src/pair_table.cpp:461)
kspace_style pppm 1.0e-4
#Neighbor style
neighbor 2.0 bin
neigh_modify check yes every 1 delay 0 page 100000 one 2000
## This part equilibrates your crystal to a pressure of ${pressure}(unit pressure) and a temperature of ${temperature}(unit temperatureture) with quantum nuclear effects
variable p_damp equal ${delta_t}*1000 #Recommended pressure damping parameter in fix nph
variable p_damp equal 0.001*1000
fix scapegoat_qtb all nph iso ${pressure} ${pressure} ${p_damp} ptemp ${temperature} #NPH does the time integration
fix scapegoat_qtb all nph iso 1.03125 ${pressure} ${p_damp} ptemp ${temperature}
fix scapegoat_qtb all nph iso 1.03125 1.03125 ${p_damp} ptemp ${temperature}
fix scapegoat_qtb all nph iso 1.03125 1.03125 1 ptemp ${temperature}
fix scapegoat_qtb all nph iso 1.03125 1.03125 1 ptemp 300
fix quartz_qtb all qtb temp ${temperature} damp ${damp_qtb} seed 35082 f_max 120.00 N_f 100 #Change f_max (THz) if your Debye frequency is higher
fix quartz_qtb all qtb temp 300 damp ${damp_qtb} seed 35082 f_max 120.00 N_f 100
fix quartz_qtb all qtb temp 300 damp 1 seed 35082 f_max 120.00 N_f 100
thermo_style custom step temp press etotal vol lx ly lz pxx pyy pzz pxy pyz pxz
thermo 200
run 2000 # 2 ps
PPPM initialization ...
using 12-bit tables for long-range coulomb (src/kspace.cpp:339)
G vector (1/distance) = 0.30159814
grid = 9 8 15
stencil order = 5
estimated absolute RMS force accuracy = 0.0011705589
estimated relative force accuracy = 8.1290814e-05
using double precision FFTW3
3d grid and FFT values/proc = 5280 1080
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 12
ghost atom cutoff = 12
binsize = 6, bins = 2 2 4
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair buck/coul/long, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
(2) pair table, perpetual, skip from (1)
attributes: half, newton on
pair build: skip
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 80.09 | 80.09 | 80.09 Mbytes
Step Temp Press TotEng Volume Lx Ly Lz Pxx Pyy Pzz Pxy Pyz Pxz
0 0 -34026.791 -2793.6042 1810.0985 9.832 8.514762 21.6216 -37470.578 -37470.432 -27139.363 3.7975984e-11 0.94245783 8.5085457e-11
200 170.7381 43248.332 -2790.8398 1879.164 9.9554912 8.6217086 21.89317 39337.624 42979.126 47428.246 324.91326 454.85872 -2034.6053
400 258.09921 -28257.8 -2788.3487 1856.1432 9.9146707 8.5863569 21.803402 -19478.873 -29571.375 -35723.151 4633.9025 8487.8103 -626.12008
600 277.77032 -22751.351 -2786.2715 1866.9783 9.9339253 8.6030319 21.845744 -21727.333 -29200.028 -17326.691 -4327.8577 -8218.4994 252.30614
800 349.8665 30508.004 -2784.2204 1873.4953 9.9454706 8.6130304 21.871134 29929.053 33562.675 28032.284 -3188.5636 12329.485 7558.5604
1000 373.67652 -18839.562 -2783.2178 1855.5937 9.9136922 8.5855095 21.80125 -18063.481 -22288.32 -16166.887 -416.09489 -10368.975 9030.4151
1200 423.34739 6846.9842 -2781.9271 1896.2131 9.9855083 8.6477041 21.959181 2147.3919 11765.847 6627.7141 -7627.6762 -1297.649 -4758.4757
1400 418.54526 -6416.7547 -2781.4358 1834.2719 9.8755745 8.5524986 21.717425 5693.0508 -19487.901 -5455.4139 827.66188 -523.1469 -3890.9904
1600 429.42798 3939.889 -2780.5861 1895.8859 9.984934 8.6472068 21.957918 3755.6972 -1326.4252 9390.395 1948.1084 4489.8536 1466.083
1800 447.76215 -8344.6447 -2780.1071 1858.4925 9.9188519 8.5899779 21.812596 -17549.502 3336.8092 -10821.241 1643.4315 -644.54621 -8935.98
2000 438.1305 -6691.4324 -2780.7407 1871.3547 9.9416812 8.6097487 21.8628 -6959.1834 -8486.8262 -4628.2877 -1019.8998 443.04638 -2751.9173
Loop time of 11.2763 on 1 procs for 2000 steps with 144 atoms
Performance: 15.324 ns/day, 1.566 hours/ns, 177.363 timesteps/s
99.5% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 7.9085 | 7.9085 | 7.9085 | 0.0 | 70.13
Kspace | 2.0339 | 2.0339 | 2.0339 | 0.0 | 18.04
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.15276 | 0.15276 | 0.15276 | 0.0 | 1.35
Output | 0.00036049 | 0.00036049 | 0.00036049 | 0.0 | 0.00
Modify | 1.1706 | 1.1706 | 1.1706 | 0.0 | 10.38
Other | | 0.01023 | | | 0.09
Nlocal: 144.000 ave 144 max 144 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 3943.00 ave 3943 max 3943 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 41952.0 ave 41952 max 41952 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 41952
Ave neighs/atom = 291.33333
Neighbor list builds = 0
Dangerous builds = 0
unfix quartz_qtb
unfix scapegoat_qtb
## This part demonstrates how to retart fix qbmsst during any stage of the shock simulation.
## PPPM may break down when compression ratio becomes extremely large. One can always use this restart technique to resume the shock simulation.
#Compression restart 1
reset_timestep 0
#Beta is the number of time steps between each update of the quantum bath temperature. Setting a larger beta can reduce thermal flactuations.
fix shock all qbmsst z ${v_msst} q ${q_msst} tscale ${tscale_msst} damp ${damp_qtb} f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q ${q_msst} tscale ${tscale_msst} damp ${damp_qtb} f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q 40 tscale ${tscale_msst} damp ${damp_qtb} f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q 40 tscale 0.05 damp ${damp_qtb} f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q 40 tscale 0.05 damp 1 f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q 40 tscale 0.05 damp 1 f_max 120 N_f 100 seed 35082 eta 1 beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q 40 tscale 0.05 damp 1 f_max 120 N_f 100 seed 35082 eta 1 beta 5 T_init 300
QBMSST parameters:
Shock in z direction
Cell mass-like parameter qmass (units of mass^2/length^4) = 4.00000e+01
Shock velocity = 7.80000e+01
Artificial viscosity (units of mass/length/time) = 0.00000e+00
Initial pressure calculated on first step
Initial volume calculated on first step
Initial energy calculated on first step
variable dhug equal f_shock[1]
variable dray equal f_shock[2]
variable lgr_vel equal f_shock[3]
variable lgr_pos equal f_shock[4]
variable T_qm equal f_shock[5] #Temperature with quantum nuclear correction
thermo_style custom step v_T_qm press econserve vol lx ly lz pzz v_dhug v_dray v_lgr_vel v_lgr_pos
thermo 200
timestep ${delta_t}
timestep 0.001
run 1000
PPPM initialization ...
using 12-bit tables for long-range coulomb (src/kspace.cpp:339)
G vector (1/distance) = 0.30313178
grid = 9 8 16
stencil order = 5
estimated absolute RMS force accuracy = 0.0010469888
estimated relative force accuracy = 7.2709348e-05
using double precision FFTW3
3d grid and FFT values/proc = 5520 1152
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 12
ghost atom cutoff = 12
binsize = 6, bins = 2 2 4
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair buck/coul/long, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
(2) pair table, perpetual, skip from (1)
attributes: half, newton on
pair build: skip
stencil: none
bin: none
Fix QBMSST v0 = 1.87135e+03
Fix QBMSST p0 = -4.62942e+03
Fix QBMSST e0 = to be -2.78074e+03
Fix QBMSST initial strain rate of -4.01095e-01 established by reducing temperature by factor of 5.00000e-02
Per MPI rank memory allocation (min/avg/max) = 80.10 | 80.10 | 80.10 Mbytes
Step v_T_qm Press Econserve Volume Lx Ly Lz Pzz v_dhug v_dray v_lgr_vel v_lgr_pos
0 300 -6922.9066 -2780.7394 1871.3547 9.9416812 8.6097487 21.8628 -4819.9374 10.953262 -190.5127 0 0
200 294.95802 54876.628 -2779.2988 1723.7617 9.9416812 8.6097487 20.13849 108897.62 -29.773363 -9271.7016 6.1518278 -15.057866
400 288.37122 139520.66 -2778.7321 1628.5573 9.9416812 8.6097487 19.02623 222107.14 8.0673735 24726.892 10.120044 -28.714689
600 280.56538 98072.818 -2779.8934 1687.2396 9.9416812 8.6097487 19.711808 164562.57 2.6099747 16006.563 7.6741039 -42.704989
800 274.9472 106058.35 -2779.2916 1651.0755 9.9416812 8.6097487 19.289307 176839.13 -39.647552 -1805.8176 9.1814643 -56.628046
1000 268.4714 189679.65 -2779.4952 1492.6558 9.9416812 8.6097487 17.43851 277332.66 -84.846841 -33118.917 15.784559 -69.870561
Loop time of 8.7779 on 1 procs for 1000 steps with 144 atoms
Performance: 9.843 ns/day, 2.438 hours/ns, 113.922 timesteps/s
99.2% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 6.8031 | 6.8031 | 6.8031 | 0.0 | 77.50
Kspace | 1.0505 | 1.0505 | 1.0505 | 0.0 | 11.97
Neigh | 0.024976 | 0.024976 | 0.024976 | 0.0 | 0.28
Comm | 0.082612 | 0.082612 | 0.082612 | 0.0 | 0.94
Output | 0.00032592 | 0.00032592 | 0.00032592 | 0.0 | 0.00
Modify | 0.8108 | 0.8108 | 0.8108 | 0.0 | 9.24
Other | | 0.005632 | | | 0.06
Nlocal: 144.000 ave 144 max 144 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 4243.00 ave 4243 max 4243 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 48210.0 ave 48210 max 48210 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 48210
Ave neighs/atom = 334.79167
Neighbor list builds = 8
Dangerous builds = 0
write_restart restart.1000
System init for write_restart ...
PPPM initialization ...
using 12-bit tables for long-range coulomb (src/kspace.cpp:339)
G vector (1/distance) = 0.30643517
grid = 9 8 15
stencil order = 5
estimated absolute RMS force accuracy = 0.0009556927
estimated relative force accuracy = 6.6369185e-05
using double precision FFTW3
3d grid and FFT values/proc = 5280 1080
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 12
ghost atom cutoff = 12
binsize = 6, bins = 2 2 3
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair buck/coul/long, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
(2) pair table, perpetual, skip from (1)
attributes: half, newton on
pair build: skip
stencil: none
bin: none
#Compression restart 2
#Read restart file and load potential again
clear
using 1 OpenMP thread(s) per MPI task
read_restart restart.1000
Reading restart file ...
restart file = 24 Dec 2020, LAMMPS = 24 Dec 2020
restoring atom style charge from restart
orthogonal box = (-0.054840605 -0.047493358 2.0915450) to (9.8868406 8.5622554 19.530055)
1 by 1 by 1 MPI processor grid
restoring pair style hybrid/overlay from restart
144 atoms
read_restart CPU = 0.001 seconds
include alpha_quartz_potential.mod
#This script implements the BKS pair potential for various silicon dioxide compounds. Inner part is fixed with a harmonic potential. Long range Coulomb interactions are evaluated with the pppm method.
#Pair Potentials
pair_style hybrid/overlay buck/coul/long ${cut_off} table linear 39901
pair_style hybrid/overlay buck/coul/long 10 table linear 39901
pair_coeff 1 1 buck/coul/long 0.0 1.0 0.0 #No interactions between Si atoms
pair_coeff 1 2 buck/coul/long 18003.757200 0.205205 133.538100
pair_coeff 2 2 buck/coul/long 1388.773000 0.362319 175.000000 #BKS interaction in PRL 64 1955 (1990)
pair_modify shift yes
pair_coeff 1 2 table potential_SiO2.TPF Si-O ${cut_off}
pair_coeff 1 2 table potential_SiO2.TPF Si-O 10
pair_coeff 2 2 table potential_SiO2.TPF O-O ${cut_off} #See the potential file for more information
pair_coeff 2 2 table potential_SiO2.TPF O-O 10
WARNING: 1 of 39901 force values in table O-O are inconsistent with -dE/dr.
Should only be flagged at inflection points (src/pair_table.cpp:461)
kspace_style pppm 1.0e-4
#Neighbor style
neighbor 2.0 bin
neigh_modify check yes every 1 delay 0 page 100000 one 2000
#Use the same fix id and add no tscale if the system is already compressed
fix shock all qbmsst z ${v_msst} q ${q_msst} tscale 0.0 damp ${damp_qtb} f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q ${q_msst} tscale 0.0 damp ${damp_qtb} f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q 40 tscale 0.0 damp ${damp_qtb} f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q 40 tscale 0.0 damp 1 f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q 40 tscale 0.0 damp 1 f_max 120 N_f 100 seed 35082 eta 1 beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q 40 tscale 0.0 damp 1 f_max 120 N_f 100 seed 35082 eta 1 beta 5 T_init 300
QBMSST parameters:
Shock in z direction
Cell mass-like parameter qmass (units of mass^2/length^4) = 4.00000e+01
Shock velocity = 7.80000e+01
Artificial viscosity (units of mass/length/time) = 0.00000e+00
Initial pressure calculated on first step
Initial volume calculated on first step
Initial energy calculated on first step
Resetting global fix info from restart file:
fix style: qbmsst, fix ID: shock
variable dhug equal f_shock[1]
variable dray equal f_shock[2]
variable lgr_vel equal f_shock[3]
variable lgr_pos equal f_shock[4]
variable T_qm equal f_shock[5] #Temperature with quantum nuclear correction
thermo_style custom step v_T_qm press econserve vol lx ly lz pzz v_dhug v_dray v_lgr_vel v_lgr_pos
thermo 500
timestep ${delta_t}
timestep 0.001
#restart 1000 restart
run 10000 #10 ps
PPPM initialization ...
using 12-bit tables for long-range coulomb (src/kspace.cpp:339)
G vector (1/distance) = 0.30643517
grid = 9 8 15
stencil order = 5
estimated absolute RMS force accuracy = 0.0009556927
estimated relative force accuracy = 6.6369185e-05
using double precision FFTW3
3d grid and FFT values/proc = 5280 1080
All restart file global fix info was re-assigned
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 12
ghost atom cutoff = 12
binsize = 6, bins = 2 2 3
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair buck/coul/long, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
(2) pair table, perpetual, skip from (1)
attributes: half, newton on
pair build: skip
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 80.12 | 80.12 | 80.12 Mbytes
Step v_T_qm Press Econserve Volume Lx Ly Lz Pzz v_dhug v_dray v_lgr_vel v_lgr_pos
1000 268.4714 189671.09 -2781.5189 1492.6558 9.9416812 8.6097487 17.43851 277359.52 -84.704915 -33092.054 15.784559 0
1500 362.24943 690690.88 -2801.3189 1013.1871 9.9416812 8.6097487 11.836938 660912.81 210.15736 -48461.278 35.76931 -24.117097
2000 851.29288 687202.75 -2829.28 998.92158 9.9416812 8.6097487 11.670275 702779.72 1177.175 -18463.457 36.363914 -45.102068
2500 1584.7231 702373.6 -2840.9926 998.47448 9.9416812 8.6097487 11.665052 685551.65 1366.7769 -36063.512 36.382549 -65.966813
3000 2369.1915 765783.04 -2835.3495 994.53954 9.9416812 8.6097487 11.61908 742616.98 1609.2532 17727.912 36.546562 -86.541922
3500 3092.3052 829161.44 -2827.8974 977.82563 9.9416812 8.6097487 11.423814 768559.05 1500.9858 29763.85 37.243215 -107.10048
4000 3627.1538 773057.81 -2813.687 988.31955 9.9416812 8.6097487 11.546413 738541.81 959.39292 8477.6525 36.805818 -127.75274
4500 3910.5455 753799.74 -2790.7863 1002.749 9.9416812 8.6097487 11.71499 763069.75 602.0127 45010.991 36.204385 -148.51594
5000 3976.7913 761978.62 -2782.0448 998.85434 9.9416812 8.6097487 11.669489 780709.6 482.34112 59410.482 36.366717 -169.42341
5500 3928.0113 702739.91 -2765.8379 1008.695 9.9416812 8.6097487 11.784456 632171.28 -991.2791 -80940.344 35.956549 -190.30951
6000 3731.5486 654300.14 -2763.1253 1032.1476 9.9416812 8.6097487 12.05845 642590.19 -832.8417 -51008.602 34.979018 -211.41573
6500 3505.5984 713092.19 -2767.3169 1010.6873 9.9416812 8.6097487 11.807732 735218.98 -176.4579 23764.995 35.873507 -232.57305
7000 3348.5047 762624.48 -2769.0996 1010.8032 9.9416812 8.6097487 11.809086 662703.98 -667.90587 -48653.562 35.868676 -253.64668
7500 3197.2839 689038.79 -2770.7583 1036.994 9.9416812 8.6097487 12.11507 679188.92 -330.37222 -10377.635 34.777016 -275.06425
8000 3117.1867 765531.79 -2775.0143 1023.9741 9.9416812 8.6097487 11.96296 681640.02 -288.82226 -18759.215 35.319699 -296.38453
8500 3053.5599 667992.24 -2772.057 1027.6458 9.9416812 8.6097487 12.005857 657921.43 -507.92809 -39422.884 35.166657 -317.64367
9000 2997.4957 704542.99 -2780.9279 1020.074 9.9416812 8.6097487 11.917396 647510.98 -398.1601 -56133.168 35.482259 -338.8598
9500 2990.5818 810181.5 -2783.2413 1002.8927 9.9416812 8.6097487 11.716669 859476.58 1107.5241 141537.42 36.198393 -359.85577
10000 3055.8298 792271.02 -2786.5277 991.62826 9.9416812 8.6097487 11.585068 847298.17 1051.5369 119986.89 36.667907 -380.83279
10500 3159.7134 706528.08 -2793.5555 1009.9173 9.9416812 8.6097487 11.798737 673106.83 67.67511 -38987.761 35.905599 -401.74903
11000 3261.7609 748345.85 -2783.0699 1017.8806 9.9416812 8.6097487 11.891771 656140.32 -267.8786 -49328.779 35.573683 -422.65593
Loop time of 130.289 on 1 procs for 10000 steps with 144 atoms
Performance: 6.631 ns/day, 3.619 hours/ns, 76.752 timesteps/s
99.6% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 109.8 | 109.8 | 109.8 | 0.0 | 84.27
Kspace | 10.328 | 10.328 | 10.328 | 0.0 | 7.93
Neigh | 1.0855 | 1.0855 | 1.0855 | 0.0 | 0.83
Comm | 1.2041 | 1.2041 | 1.2041 | 0.0 | 0.92
Output | 0.0012848 | 0.0012848 | 0.0012848 | 0.0 | 0.00
Modify | 7.8094 | 7.8094 | 7.8094 | 0.0 | 5.99
Other | | 0.06511 | | | 0.05
Nlocal: 144.000 ave 144 max 144 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 5430.00 ave 5430 max 5430 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 72807.0 ave 72807 max 72807 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 72807
Ave neighs/atom = 505.60417
Neighbor list builds = 206
Dangerous builds = 0
shell rm restart.1000
Total wall time: 0:02:32

View File

@ -0,0 +1,435 @@
LAMMPS (24 Dec 2020)
using 1 OpenMP thread(s) per MPI task
## This script first uses fix qtb to equilibrate alpha quartz structure to an initial state with quantum nuclear correction and then simulate shock induced phase transition through the quantum thermal bath multi-scale shock technique
variable x_rep equal 2 #plot is made with x_rep = 8 #x-direction replication number
variable y_rep equal 1 #plot is made with y_rep = 5 #y-direction replication number
variable z_rep equal 4 #plot is made with z_rep = 15 #z-direction replication number
variable cut_off equal 10.0 #Cut-off distance for the Buckingham term (Angstrom in metal units)
variable pressure equal 1.03125 #Initial state pressure (bar in metal units)
variable temperature equal 300.0 #Initial state quantum temperature (K in metal units)
variable delta_t equal 1.0e-3 #MD timestep length (ps in metal units)
variable damp_qtb equal 1 #1/gamma where gamma is the friction coefficient in quantum thermal bath (ps in metal units)
variable v_msst equal 78.0 #Shock velocity (Angstrom/ps in metal units)
variable q_msst equal 40.0 #Box mass-like parameter in the MSST (mass^2/length^4, where mass=grams/mole and length=Angstrom in metal units)
variable tscale_msst equal 0.05 #Temperature reduction parameter in the MSST (unitless)
variable eta_qbmsst equal 1.0 #Coupling constant between the shock and the quantum thermal bath (unitless constant)
## This part uses fix qtb to prepare alpha-quartz with quantum nuclear correction of the initial state
include alpha_quartz_qtb.mod
## This script first constructs an alpha quartz structure of a given size. It then uses fix qtb to equilibrate the computational cell to the specified temperature and pressure.
## This part defines units, alpha-quartz crystal, and atomic information
#General
units metal
dimension 3
boundary p p p
atom_style charge
#Lattice
lattice custom 1.0 a1 4.916000 0.000000 0.000000 a2 -2.45800 4.257381 0.000000 a3 0.000000 0.000000 5.405400 basis 0.469700 0.000000 0.000000 basis 0.000000 0.469700 0.666667 basis 0.530300 0.530300 0.333333 basis 0.413500 0.266900 0.119100 basis 0.266900 0.413500 0.547567 basis 0.733100 0.146600 0.785767 basis 0.586500 0.853400 0.214233 basis 0.853400 0.586500 0.452433 basis 0.146600 0.733100 0.880900 #American Mineralogist 65 920 1980 (Space Group 154)
Lattice spacing in x,y,z = 7.3740000 4.2573810 5.4054000
#Computational Cell
region orthorhombic_unit_cell block 0 4.916000 0 8.514762 0 5.405400 units box
create_box 2 orthorhombic_unit_cell
Created orthogonal box = (0.0000000 0.0000000 0.0000000) to (4.9160000 8.5147620 5.4054000)
1 by 2 by 2 MPI processor grid
create_atoms 1 box basis 1 1 basis 2 1 basis 3 1 basis 4 2 basis 5 2 basis 6 2 basis 7 2 basis 8 2 basis 9 2
Created 18 atoms
create_atoms CPU = 0.000 seconds
replicate ${x_rep} ${y_rep} ${z_rep}
replicate 2 ${y_rep} ${z_rep}
replicate 2 1 ${z_rep}
replicate 2 1 4
Replicating atoms ...
orthogonal box = (0.0000000 0.0000000 0.0000000) to (9.8320000 8.5147620 21.621600)
1 by 1 by 4 MPI processor grid
144 atoms
replicate CPU = 0.001 seconds
#Atomic Information
mass 1 28.085500
mass 2 15.999400
set type 1 charge +2.4
Setting atom values ...
48 settings made for charge
set type 2 charge -1.2
Setting atom values ...
96 settings made for charge
## This part implements the BKS pair potential with a cut-off distance for the Buckingham term. Long range Coulomb interactions are evaluated with the pppm method.
include alpha_quartz_potential.mod
#This script implements the BKS pair potential for various silicon dioxide compounds. Inner part is fixed with a harmonic potential. Long range Coulomb interactions are evaluated with the pppm method.
#Pair Potentials
pair_style hybrid/overlay buck/coul/long ${cut_off} table linear 39901
pair_style hybrid/overlay buck/coul/long 10 table linear 39901
pair_coeff 1 1 buck/coul/long 0.0 1.0 0.0 #No interactions between Si atoms
pair_coeff 1 2 buck/coul/long 18003.757200 0.205205 133.538100
pair_coeff 2 2 buck/coul/long 1388.773000 0.362319 175.000000 #BKS interaction in PRL 64 1955 (1990)
pair_modify shift yes
pair_coeff 1 2 table potential_SiO2.TPF Si-O ${cut_off}
pair_coeff 1 2 table potential_SiO2.TPF Si-O 10
pair_coeff 2 2 table potential_SiO2.TPF O-O ${cut_off} #See the potential file for more information
pair_coeff 2 2 table potential_SiO2.TPF O-O 10
WARNING: 1 of 39901 force values in table O-O are inconsistent with -dE/dr.
Should only be flagged at inflection points (src/pair_table.cpp:461)
kspace_style pppm 1.0e-4
#Neighbor style
neighbor 2.0 bin
neigh_modify check yes every 1 delay 0 page 100000 one 2000
## This part equilibrates your crystal to a pressure of ${pressure}(unit pressure) and a temperature of ${temperature}(unit temperatureture) with quantum nuclear effects
variable p_damp equal ${delta_t}*1000 #Recommended pressure damping parameter in fix nph
variable p_damp equal 0.001*1000
fix scapegoat_qtb all nph iso ${pressure} ${pressure} ${p_damp} ptemp ${temperature} #NPH does the time integration
fix scapegoat_qtb all nph iso 1.03125 ${pressure} ${p_damp} ptemp ${temperature}
fix scapegoat_qtb all nph iso 1.03125 1.03125 ${p_damp} ptemp ${temperature}
fix scapegoat_qtb all nph iso 1.03125 1.03125 1 ptemp ${temperature}
fix scapegoat_qtb all nph iso 1.03125 1.03125 1 ptemp 300
fix quartz_qtb all qtb temp ${temperature} damp ${damp_qtb} seed 35082 f_max 120.00 N_f 100 #Change f_max (THz) if your Debye frequency is higher
fix quartz_qtb all qtb temp 300 damp ${damp_qtb} seed 35082 f_max 120.00 N_f 100
fix quartz_qtb all qtb temp 300 damp 1 seed 35082 f_max 120.00 N_f 100
thermo_style custom step temp press etotal vol lx ly lz pxx pyy pzz pxy pyz pxz
thermo 200
run 2000 # 2 ps
PPPM initialization ...
using 12-bit tables for long-range coulomb (src/kspace.cpp:339)
G vector (1/distance) = 0.30159814
grid = 9 8 15
stencil order = 5
estimated absolute RMS force accuracy = 0.0011705589
estimated relative force accuracy = 8.1290814e-05
using double precision FFTW3
3d grid and FFT values/proc = 2400 288
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 12
ghost atom cutoff = 12
binsize = 6, bins = 2 2 4
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair buck/coul/long, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
(2) pair table, perpetual, skip from (1)
attributes: half, newton on
pair build: skip
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 79.70 | 79.70 | 79.70 Mbytes
Step Temp Press TotEng Volume Lx Ly Lz Pxx Pyy Pzz Pxy Pyz Pxz
0 0 -34026.791 -2793.6042 1810.0985 9.832 8.514762 21.6216 -37470.578 -37470.432 -27139.363 -3.2012455e-11 0.94245783 1.6892124e-10
200 153.57631 45538.205 -2790.8177 1873.0866 9.9447472 8.612404 21.869543 41721.016 44095.248 50798.351 -3961.4596 1223.325 2871.656
400 234.74785 -34404.175 -2789.0189 1850.2127 9.9041 8.5772024 21.780156 -28329.333 -39376.313 -35506.88 -1154.5043 -5411.1071 2246.6749
600 265.24834 -20905.145 -2786.2727 1874.9981 9.948129 8.6153326 21.87698 -22753.885 -21091.083 -18870.467 -4645.5539 2968.2936 1415.0335
800 297.79036 32990.577 -2784.8247 1853.6946 9.910309 8.5825796 21.79381 30061.366 35359.175 33551.191 -3092.2938 1525.518 -6461.029
1000 367.71885 -27539.237 -2783.0102 1864.7161 9.9299114 8.5995557 21.836917 -20273.384 -38720.43 -23623.895 7639.0325 -866.34777 543.5312
1200 399.7711 3807.785 -2781.511 1893.4978 9.9807399 8.6435745 21.948695 1625.8297 7441.2317 2356.2937 -4057.1659 3814.9292 1528.4637
1400 466.57958 -4148.2231 -2780.1546 1851.5925 9.9065614 8.5793341 21.785568 -10883.182 1816.778 -3378.2653 896.24645 -7208.5417 -42.262464
1600 497.86536 14505.308 -2778.9409 1882.2616 9.9609584 8.6264432 21.905193 8268.1088 20614.74 14633.075 -2690.5703 6807.3188 11995.875
1800 557.31178 -108.02787 -2778.1875 1875.514 9.9490413 8.6161228 21.878986 948.70277 -1929.753 656.96663 -1628.2124 -6594.6026 -4423.4256
2000 480.39444 -8852.2282 -2778.4963 1862.9552 9.9267847 8.596848 21.830042 -18274.302 3038.8276 -11321.21 -5002.1095 12023.298 6845.2631
Loop time of 4.1373 on 4 procs for 2000 steps with 144 atoms
Performance: 41.766 ns/day, 0.575 hours/ns, 483.407 timesteps/s
96.8% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 2.2821 | 2.4503 | 2.7881 | 12.7 | 59.22
Kspace | 0.81032 | 1.1413 | 1.306 | 18.2 | 27.59
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.16143 | 0.16964 | 0.17659 | 1.3 | 4.10
Output | 0.00026584 | 0.00061899 | 0.0016773 | 0.0 | 0.01
Modify | 0.29943 | 0.33639 | 0.34927 | 3.7 | 8.13
Other | | 0.03911 | | | 0.95
Nlocal: 36.0000 ave 36 max 36 min
Histogram: 4 0 0 0 0 0 0 0 0 0
Nghost: 2614.00 ave 2614 max 2614 min
Histogram: 4 0 0 0 0 0 0 0 0 0
Neighs: 10488.0 ave 11326 max 9404 min
Histogram: 1 0 0 0 0 0 2 0 0 1
Total # of neighbors = 41952
Ave neighs/atom = 291.33333
Neighbor list builds = 0
Dangerous builds = 0
unfix quartz_qtb
unfix scapegoat_qtb
## This part demonstrates how to retart fix qbmsst during any stage of the shock simulation.
## PPPM may break down when compression ratio becomes extremely large. One can always use this restart technique to resume the shock simulation.
#Compression restart 1
reset_timestep 0
#Beta is the number of time steps between each update of the quantum bath temperature. Setting a larger beta can reduce thermal flactuations.
fix shock all qbmsst z ${v_msst} q ${q_msst} tscale ${tscale_msst} damp ${damp_qtb} f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q ${q_msst} tscale ${tscale_msst} damp ${damp_qtb} f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q 40 tscale ${tscale_msst} damp ${damp_qtb} f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q 40 tscale 0.05 damp ${damp_qtb} f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q 40 tscale 0.05 damp 1 f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q 40 tscale 0.05 damp 1 f_max 120 N_f 100 seed 35082 eta 1 beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q 40 tscale 0.05 damp 1 f_max 120 N_f 100 seed 35082 eta 1 beta 5 T_init 300
QBMSST parameters:
Shock in z direction
Cell mass-like parameter qmass (units of mass^2/length^4) = 4.00000e+01
Shock velocity = 7.80000e+01
Artificial viscosity (units of mass/length/time) = 0.00000e+00
Initial pressure calculated on first step
Initial volume calculated on first step
Initial energy calculated on first step
variable dhug equal f_shock[1]
variable dray equal f_shock[2]
variable lgr_vel equal f_shock[3]
variable lgr_pos equal f_shock[4]
variable T_qm equal f_shock[5] #Temperature with quantum nuclear correction
thermo_style custom step v_T_qm press econserve vol lx ly lz pzz v_dhug v_dray v_lgr_vel v_lgr_pos
thermo 200
timestep ${delta_t}
timestep 0.001
run 1000
PPPM initialization ...
using 12-bit tables for long-range coulomb (src/kspace.cpp:339)
G vector (1/distance) = 0.30087967
grid = 9 8 15
stencil order = 5
estimated absolute RMS force accuracy = 0.0012053392
estimated relative force accuracy = 8.3706174e-05
using double precision FFTW3
3d grid and FFT values/proc = 2400 288
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 12
ghost atom cutoff = 12
binsize = 6, bins = 2 2 4
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair buck/coul/long, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
(2) pair table, perpetual, skip from (1)
attributes: half, newton on
pair build: skip
stencil: none
bin: none
Fix QBMSST v0 = 1.86296e+03
Fix QBMSST p0 = -1.13219e+04
Fix QBMSST e0 = to be -2.77850e+03
Fix QBMSST initial strain rate of -4.21889e-01 established by reducing temperature by factor of 5.00000e-02
Per MPI rank memory allocation (min/avg/max) = 79.70 | 79.70 | 79.70 Mbytes
Step v_T_qm Press Econserve Volume Lx Ly Lz Pzz v_dhug v_dray v_lgr_vel v_lgr_pos
0 300 -9106.3219 -2778.4963 1862.9552 9.9267847 8.596848 21.830042 -11562.009 12.009861 -240.06987 0 0
200 296.47212 25984.099 -2777.5178 1770.2165 9.9267847 8.596848 20.743332 64970.178 -25.305804 -1564.7427 3.8828751 -15.16768
400 291.06704 69977.415 -2777.6325 1684.8932 9.9267847 8.596848 19.743517 144833.61 -12.18477 6667.3264 7.4552723 -29.607029
600 287.21114 39706.769 -2778.0322 1716.9533 9.9267847 8.596848 20.119196 87971.211 -38.594057 -23279.705 6.1129499 -43.7513
800 284.33606 18833.325 -2778.1637 1792.7575 9.9267847 8.596848 21.007467 43725.516 -8.1270751 -3885.5508 2.9391052 -58.454557
1000 281.98323 -6030.7047 -2778.3314 1881.8368 9.9267847 8.596848 22.051295 -14118.589 1.3182589 13054.989 -0.79055248 -73.780966
Loop time of 3.32539 on 4 procs for 1000 steps with 144 atoms
Performance: 25.982 ns/day, 0.924 hours/ns, 300.717 timesteps/s
97.0% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 1.9626 | 2.0842 | 2.2541 | 7.9 | 62.68
Kspace | 0.44255 | 0.61231 | 0.73369 | 14.5 | 18.41
Neigh | 0.0050733 | 0.0052404 | 0.0053804 | 0.2 | 0.16
Comm | 0.077084 | 0.077385 | 0.077714 | 0.1 | 2.33
Output | 0.00029039 | 0.00046909 | 0.0010037 | 0.0 | 0.01
Modify | 0.50853 | 0.52962 | 0.53724 | 1.7 | 15.93
Other | | 0.01615 | | | 0.49
Nlocal: 36.0000 ave 38 max 34 min
Histogram: 1 0 1 0 0 0 0 1 0 1
Nghost: 2527.75 ave 2547 max 2518 min
Histogram: 2 0 0 1 0 0 0 0 0 1
Neighs: 10194.8 ave 11177 max 9437 min
Histogram: 2 0 0 0 0 0 1 0 0 1
Total # of neighbors = 40779
Ave neighs/atom = 283.18750
Neighbor list builds = 6
Dangerous builds = 0
write_restart restart.1000
System init for write_restart ...
PPPM initialization ...
using 12-bit tables for long-range coulomb (src/kspace.cpp:339)
G vector (1/distance) = 0.30295266
grid = 9 8 16
stencil order = 5
estimated absolute RMS force accuracy = 0.0010556863
estimated relative force accuracy = 7.3313358e-05
using double precision FFTW3
3d grid and FFT values/proc = 2640 288
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 12
ghost atom cutoff = 12
binsize = 6, bins = 2 2 4
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair buck/coul/long, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
(2) pair table, perpetual, skip from (1)
attributes: half, newton on
pair build: skip
stencil: none
bin: none
#Compression restart 2
#Read restart file and load potential again
clear
using 1 OpenMP thread(s) per MPI task
read_restart restart.1000
Reading restart file ...
restart file = 24 Dec 2020, LAMMPS = 24 Dec 2020
restoring atom style charge from restart
orthogonal box = (-0.047392358 -0.041042987 -0.21484765) to (9.8793924 8.5558050 21.836448)
1 by 1 by 4 MPI processor grid
restoring pair style hybrid/overlay from restart
144 atoms
read_restart CPU = 0.009 seconds
include alpha_quartz_potential.mod
#This script implements the BKS pair potential for various silicon dioxide compounds. Inner part is fixed with a harmonic potential. Long range Coulomb interactions are evaluated with the pppm method.
#Pair Potentials
pair_style hybrid/overlay buck/coul/long ${cut_off} table linear 39901
pair_style hybrid/overlay buck/coul/long 10 table linear 39901
pair_coeff 1 1 buck/coul/long 0.0 1.0 0.0 #No interactions between Si atoms
pair_coeff 1 2 buck/coul/long 18003.757200 0.205205 133.538100
pair_coeff 2 2 buck/coul/long 1388.773000 0.362319 175.000000 #BKS interaction in PRL 64 1955 (1990)
pair_modify shift yes
pair_coeff 1 2 table potential_SiO2.TPF Si-O ${cut_off}
pair_coeff 1 2 table potential_SiO2.TPF Si-O 10
pair_coeff 2 2 table potential_SiO2.TPF O-O ${cut_off} #See the potential file for more information
pair_coeff 2 2 table potential_SiO2.TPF O-O 10
WARNING: 1 of 39901 force values in table O-O are inconsistent with -dE/dr.
Should only be flagged at inflection points (src/pair_table.cpp:461)
kspace_style pppm 1.0e-4
#Neighbor style
neighbor 2.0 bin
neigh_modify check yes every 1 delay 0 page 100000 one 2000
#Use the same fix id and add no tscale if the system is already compressed
fix shock all qbmsst z ${v_msst} q ${q_msst} tscale 0.0 damp ${damp_qtb} f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q ${q_msst} tscale 0.0 damp ${damp_qtb} f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q 40 tscale 0.0 damp ${damp_qtb} f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q 40 tscale 0.0 damp 1 f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q 40 tscale 0.0 damp 1 f_max 120 N_f 100 seed 35082 eta 1 beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q 40 tscale 0.0 damp 1 f_max 120 N_f 100 seed 35082 eta 1 beta 5 T_init 300
QBMSST parameters:
Shock in z direction
Cell mass-like parameter qmass (units of mass^2/length^4) = 4.00000e+01
Shock velocity = 7.80000e+01
Artificial viscosity (units of mass/length/time) = 0.00000e+00
Initial pressure calculated on first step
Initial volume calculated on first step
Initial energy calculated on first step
Resetting global fix info from restart file:
fix style: qbmsst, fix ID: shock
variable dhug equal f_shock[1]
variable dray equal f_shock[2]
variable lgr_vel equal f_shock[3]
variable lgr_pos equal f_shock[4]
variable T_qm equal f_shock[5] #Temperature with quantum nuclear correction
thermo_style custom step v_T_qm press econserve vol lx ly lz pzz v_dhug v_dray v_lgr_vel v_lgr_pos
thermo 500
timestep ${delta_t}
timestep 0.001
#restart 1000 restart
run 10000 #10 ps
PPPM initialization ...
using 12-bit tables for long-range coulomb (src/kspace.cpp:339)
G vector (1/distance) = 0.30295266
grid = 9 8 16
stencil order = 5
estimated absolute RMS force accuracy = 0.0010556863
estimated relative force accuracy = 7.3313358e-05
using double precision FFTW3
3d grid and FFT values/proc = 2640 288
All restart file global fix info was re-assigned
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 12
ghost atom cutoff = 12
binsize = 6, bins = 2 2 4
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair buck/coul/long, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
(2) pair table, perpetual, skip from (1)
attributes: half, newton on
pair build: skip
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 79.71 | 79.71 | 79.71 Mbytes
Step v_T_qm Press Econserve Volume Lx Ly Lz Pzz v_dhug v_dray v_lgr_vel v_lgr_pos
1000 281.98323 -6031.2507 -2778.6227 1881.8368 9.9267847 8.596848 22.051295 -14113.608 1.3371988 13059.97 -0.79055248 0
1500 266.12743 44405.252 -2777.9815 1739.6551 9.9267847 8.596848 20.385215 92589.619 -12.060756 397.55607 5.1624473 -37.823753
2000 255.79412 17620.89 -2777.9685 1785.7605 9.9267847 8.596848 20.925477 48671.42 -16.082485 -4813.8454 3.2320631 -73.974438
2500 257.13592 39692.462 -2778.6986 1751.4095 9.9267847 8.596848 20.522952 80667.315 15.746345 -1656.6275 4.6703047 -112.35088
3000 248.95332 9617.5633 -2778.937 1830.5557 9.9267847 8.596848 21.450385 25275.769 19.730704 9397.3972 1.3565331 -148.37113
3500 247.70025 100159.87 -2778.0604 1610.8047 9.9267847 8.596848 18.875351 189849.69 -33.726976 -10516.027 10.557281 -185.61862
4000 266.07224 848367.31 -2787.9052 992.46097 9.9267847 8.596848 11.629622 880163.37 1477.3994 160680.23 36.44669 -213.83067
4500 645.86948 789169.63 -2822.9559 992.40405 9.9267847 8.596848 11.628955 696879.41 1039.4139 -22651.518 36.449073 -234.79958
5000 1369.4257 735014.89 -2838.4571 1002.6048 9.9267847 8.596848 11.748487 648785.76 1170.3517 -62181.314 36.021977 -255.55776
5500 2156.7632 768865.28 -2835.9297 995.94989 9.9267847 8.596848 11.670505 678013.94 1271.734 -38540.152 36.300612 -276.42588
6000 2864.2837 773631.53 -2828.0627 993.01727 9.9267847 8.596848 11.63614 749067.81 1567.7659 30051.708 36.423398 -297.26898
6500 3422.632 861319.73 -2810.1415 985.48363 9.9267847 8.596848 11.547861 816792.18 1535.8348 91451.363 36.738824 -318.12934
7000 3798.2073 791521.73 -2801.7757 993.1961 9.9267847 8.596848 11.638236 677215.78 330.09854 -41650.204 36.415911 -338.86015
7500 4060.7728 836165.25 -2789.6215 984.13658 9.9267847 8.596848 11.532077 780101.5 698.84908 53629.791 36.795223 -359.64284
8000 4122.5641 754871.86 -2776.0049 1006.6266 9.9267847 8.596848 11.795613 699610.84 -124.86381 -7979.8848 35.853592 -380.58907
8500 4087.3529 769727.63 -2775.3629 1018.2197 9.9267847 8.596848 11.931461 767853.09 415.9984 69995.141 35.368199 -401.90058
9000 3958.4459 615996.33 -2758.7864 1058.0696 9.9267847 8.596848 12.398422 641295.34 -689.82578 -23107.426 33.699723 -423.43203
9500 3746.2013 643366.31 -2767.1851 1043.1232 9.9267847 8.596848 12.22328 610176.19 -767.67823 -66774.534 34.325515 -445.14544
10000 3723.8623 659730.11 -2781.6634 1034.0441 9.9267847 8.596848 12.116891 671355.25 0.037615796 -13217.642 34.705647 -466.9448
10500 3705.48 637406.18 -2776.4898 1041.5851 9.9267847 8.596848 12.205256 725619.7 274.78304 47377.665 34.389914 -488.75102
11000 3678.0139 648116.35 -2779.0968 1049.9523 9.9267847 8.596848 12.303303 723144.21 382.51198 51926.71 34.039587 -510.63944
Loop time of 51.3151 on 4 procs for 10000 steps with 144 atoms
Performance: 16.837 ns/day, 1.425 hours/ns, 194.874 timesteps/s
94.1% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 26.662 | 32.362 | 38.05 | 70.9 | 63.07
Kspace | 5.8733 | 11.582 | 17.302 | 118.9 | 22.57
Neigh | 0.18541 | 0.22229 | 0.25113 | 5.0 | 0.43
Comm | 1.4273 | 1.4501 | 1.483 | 1.9 | 2.83
Output | 0.0011935 | 0.0018681 | 0.003891 | 2.7 | 0.00
Modify | 5.4539 | 5.5056 | 5.5294 | 1.3 | 10.73
Other | | 0.1916 | | | 0.37
Nlocal: 36.0000 ave 37 max 35 min
Histogram: 1 0 0 0 0 2 0 0 0 1
Nghost: 4159.50 ave 4171 max 4140 min
Histogram: 1 0 0 0 0 0 1 0 1 1
Neighs: 17967.8 ave 20291 max 15710 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Total # of neighbors = 71871
Ave neighs/atom = 499.10417
Neighbor list builds = 161
Dangerous builds = 0
shell rm restart.1000
Total wall time: 0:01:01

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,79 @@
## This script first constructs an alpha quartz structure of a given size. It then uses fix qtb to equilibrate the computational cell to the specified temperature and pressure.
variable x_rep equal 2 #x-direction replication number
variable y_rep equal 2 #y-direction replication number
variable z_rep equal 2 #z-direction replication number
variable cut_off equal 10.0 #Cut-off distance for the Buckingham term (Angstrom in metal units)
variable pressure equal 1.03125 #Initial state pressure (bar in metal units)
variable temperature equal 300.0 #Initial state quantum temperature (K in metal units)
variable delta_t equal 1.0e-3 #MD timestep length (ps in metal units)
variable damp_qtb equal 1 #1/gamma where gamma is the friction coefficient in quantum thermal bath (ps in metal units)
## This part defines units, alpha-quartz crystal, and atomic information
#General
units metal
dimension 3
boundary p p p
atom_style charge
#Lattice
lattice custom 1.0 &
a1 4.916000 0.000000 0.000000 &
a2 -2.45800 4.257381 0.000000 &
a3 0.000000 0.000000 5.405400 &
&
basis 0.469700 0.000000 0.000000 &
basis 0.000000 0.469700 0.666667 &
basis 0.530300 0.530300 0.333333 &
&
basis 0.413500 0.266900 0.119100 &
basis 0.266900 0.413500 0.547567 &
basis 0.733100 0.146600 0.785767 &
basis 0.586500 0.853400 0.214233 &
basis 0.853400 0.586500 0.452433 &
basis 0.146600 0.733100 0.880900 #American Mineralogist 65 920 1980 (Space Group 154)
#Computational Cell
region orthorhombic_unit_cell block 0 4.916000 0 8.514762 0 5.405400 units box
create_box 2 orthorhombic_unit_cell
create_atoms 1 box &
basis 1 1 &
basis 2 1 &
basis 3 1 &
basis 4 2 &
basis 5 2 &
basis 6 2 &
basis 7 2 &
basis 8 2 &
basis 9 2
replicate ${x_rep} ${y_rep} ${z_rep}
#Atomic Information
mass 1 28.085500
mass 2 15.999400
set type 1 charge +2.4
set type 2 charge -1.2
## This part implements the BKS pair potential with a cut-off distance for the Buckingham term. Long range Coulomb interactions are evaluated with the pppm method.
#Pair Potentials
pair_style buck/coul/long ${cut_off} #BKS interaction, PRL 64 1955 (1990)
pair_coeff 1 1 0.0 1.0 0.0
pair_coeff 1 2 18003.757200 0.205205 133.538100
pair_coeff 2 2 1388.773000 0.362319 175.000000
pair_modify shift yes
kspace_style pppm 1.0e-4
#Neighbor style
neighbor 2.0 bin
neigh_modify check yes every 1 delay 0 page 100000 one 2000
## This part equilibrates your crystal to a pressure of ${pressure}(unit pressure) and a temperature of ${temperature}(unit temperatureture) with quantum nuclear effects
variable p_damp equal ${delta_t}*1000 #Recommended pressure damping parameter in fix nph
fix scapegoat_qtb all nph iso ${pressure} ${pressure} ${p_damp} ptemp ${temperature} #NPH does the time integration
fix quartz_qtb all qtb temp ${temperature} damp ${damp_qtb} seed 35082 f_max 120.00 N_f 100 #Change f_max (THz) if your Debye frequency is higher
thermo_style custom step temp press etotal vol lx ly lz pxx pyy pzz pxy pyz pxz
thermo 500
run 10000 # 20 ps
unfix quartz_qtb
unfix scapegoat_qtb

View File

@ -0,0 +1,152 @@
LAMMPS (15 Jun 2020)
using 1 OpenMP thread(s) per MPI task
## This script first constructs an alpha quartz structure of a given size. It then uses fix qtb to equilibrate the computational cell to the specified temperature and pressure.
variable x_rep equal 2 #x-direction replication number
variable y_rep equal 2 #y-direction replication number
variable z_rep equal 2 #z-direction replication number
variable cut_off equal 10.0 #Cut-off distance for the Buckingham term (Angstrom in metal units)
variable pressure equal 1.03125 #Initial state pressure (bar in metal units)
variable temperature equal 300.0 #Initial state quantum temperature (K in metal units)
variable delta_t equal 1.0e-3 #MD timestep length (ps in metal units)
variable damp_qtb equal 1 #1/gamma where gamma is the friction coefficient in quantum thermal bath (ps in metal units)
## This part defines units, alpha-quartz crystal, and atomic information
#General
units metal
dimension 3
boundary p p p
atom_style charge
#Lattice
lattice custom 1.0 a1 4.916000 0.000000 0.000000 a2 -2.45800 4.257381 0.000000 a3 0.000000 0.000000 5.405400 basis 0.469700 0.000000 0.000000 basis 0.000000 0.469700 0.666667 basis 0.530300 0.530300 0.333333 basis 0.413500 0.266900 0.119100 basis 0.266900 0.413500 0.547567 basis 0.733100 0.146600 0.785767 basis 0.586500 0.853400 0.214233 basis 0.853400 0.586500 0.452433 basis 0.146600 0.733100 0.880900 #American Mineralogist 65 920 1980 (Space Group 154)
Lattice spacing in x,y,z = 7.374 4.25738 5.4054
#Computational Cell
region orthorhombic_unit_cell block 0 4.916000 0 8.514762 0 5.405400 units box
create_box 2 orthorhombic_unit_cell
Created orthogonal box = (0.0 0.0 0.0) to (4.916 8.514762 5.4054)
1 by 1 by 1 MPI processor grid
create_atoms 1 box basis 1 1 basis 2 1 basis 3 1 basis 4 2 basis 5 2 basis 6 2 basis 7 2 basis 8 2 basis 9 2
Created 18 atoms
create_atoms CPU = 0.000 seconds
replicate ${x_rep} ${y_rep} ${z_rep}
replicate 2 ${y_rep} ${z_rep}
replicate 2 2 ${z_rep}
replicate 2 2 2
orthogonal box = (0.0 0.0 0.0) to (9.832 17.029524 10.8108)
1 by 1 by 1 MPI processor grid
144 atoms
replicate CPU = 0.000219584 secs
#Atomic Information
mass 1 28.085500
mass 2 15.999400
set type 1 charge +2.4
48 settings made for charge
set type 2 charge -1.2
96 settings made for charge
## This part implements the BKS pair potential with a cut-off distance for the Buckingham term. Long range Coulomb interactions are evaluated with the pppm method.
#Pair Potentials
pair_style buck/coul/long ${cut_off} #BKS interaction, PRL 64 1955 (1990)
pair_style buck/coul/long 10
pair_coeff 1 1 0.0 1.0 0.0
pair_coeff 1 2 18003.757200 0.205205 133.538100
pair_coeff 2 2 1388.773000 0.362319 175.000000
pair_modify shift yes
kspace_style pppm 1.0e-4
#Neighbor style
neighbor 2.0 bin
neigh_modify check yes every 1 delay 0 page 100000 one 2000
## This part equilibrates your crystal to a pressure of ${pressure}(unit pressure) and a temperature of ${temperature}(unit temperatureture) with quantum nuclear effects
variable p_damp equal ${delta_t}*1000 #Recommended pressure damping parameter in fix nph
variable p_damp equal 0.001*1000
fix scapegoat_qtb all nph iso ${pressure} ${pressure} ${p_damp} ptemp ${temperature} #NPH does the time integration
fix scapegoat_qtb all nph iso 1.03125 ${pressure} ${p_damp} ptemp ${temperature}
fix scapegoat_qtb all nph iso 1.03125 1.03125 ${p_damp} ptemp ${temperature}
fix scapegoat_qtb all nph iso 1.03125 1.03125 1 ptemp ${temperature}
fix scapegoat_qtb all nph iso 1.03125 1.03125 1 ptemp 300
fix quartz_qtb all qtb temp ${temperature} damp ${damp_qtb} seed 35082 f_max 120.00 N_f 100 #Change f_max (THz) if your Debye frequency is higher
fix quartz_qtb all qtb temp 300 damp ${damp_qtb} seed 35082 f_max 120.00 N_f 100
fix quartz_qtb all qtb temp 300 damp 1 seed 35082 f_max 120.00 N_f 100
thermo_style custom step temp press etotal vol lx ly lz pxx pyy pzz pxy pyz pxz
thermo 500
run 10000 # 20 ps
PPPM initialization ...
using 12-bit tables for long-range coulomb (src/kspace.cpp:332)
G vector (1/distance) = 0.307414
grid = 9 15 10
stencil order = 5
estimated absolute RMS force accuracy = 0.000822922
estimated relative force accuracy = 5.71487e-05
using double precision FFTW3
3d grid and FFT values/proc = 5984 1350
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 12
ghost atom cutoff = 12
binsize = 6, bins = 2 3 2
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair buck/coul/long, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
Per MPI rank memory allocation (min/avg/max) = 79.54 | 79.54 | 79.54 Mbytes
Step Temp Press TotEng Volume Lx Ly Lz Pxx Pyy Pzz Pxy Pyz Pxz
0 0 -34025.794 -2793.6041 1810.0985 9.832 17.029524 10.8108 -37478.502 -37477.413 -27121.466 -1.3649088e-10 1.3388978 5.8209479e-10
500 281.29079 -40385.348 -2786.6755 1844.5575 9.893999 17.136909 10.878971 -44649.574 -45631.516 -30874.953 -5970.3691 3630.1324 5208.8966
1000 405.39848 -15491.657 -2783.3315 1874.0851 9.9465141 17.227868 10.936714 -7770.4561 -21469.887 -17234.627 -4706.5632 -8313.9522 -5109.7918
1500 428.98568 -12118.951 -2781.3232 1874.1627 9.9466513 17.228106 10.936865 -11239.135 -11740.052 -13377.666 -3778.9317 -6220.1431 12775.412
2000 522.11905 -6687.482 -2779.2181 1855.9626 9.914349 17.172157 10.901347 -8016.0133 -10737.23 -1309.2028 -4980.3805 5270.2674 5848.5479
2500 496.74376 4264.2623 -2778.9979 1849.3244 9.9025147 17.151659 10.888335 -477.1374 3487.19 9782.7343 -4961.2016 2380.6522 4736.0758
3000 456.49628 2320.781 -2779.3844 1853.2925 9.9095923 17.163918 10.896117 5479.6232 -2954.3023 4437.022 3719.9287 4445.0723 -3278.5058
3500 485.20722 -7480.1789 -2778.6062 1859.6305 9.920876 17.183462 10.908524 -9340.6334 -6129.8494 -6970.0541 -7379.3507 1772.8159 334.33057
4000 527.61216 -13499.73 -2777.3433 1889.9405 9.9744857 17.276316 10.96747 -16483.038 -7465.2297 -16550.923 -2517.02 -1863.063 3314.927
4500 519.94117 721.60614 -2777.8506 1879.6562 9.9563603 17.244922 10.947541 -913.2791 -1765.7541 4843.8516 4466.5704 -14141.087 -6439.5669
5000 505.27757 -6278.3805 -2777.3641 1881.2931 9.9592497 17.249927 10.950718 -14254.233 -2653.6233 -1927.2858 1838.1568 5767.9267 597.47761
5500 500.70903 11303.821 -2777.8881 1871.0076 9.9410666 17.218433 10.930724 -6452.7947 24876.967 15487.29 522.01171 10473.257 9780.893
6000 526.65329 7991.2419 -2777.172 1856.9227 9.9160583 17.175117 10.903227 -68.823156 11005.468 13037.081 1253.9214 10039.559 1053.0486
6500 485.30026 12811.546 -2777.5866 1845.31 9.8953442 17.139239 10.88045 10063.921 20215.037 8155.6798 -3886.954 2249.2807 4855.0011
7000 507.85472 2649.7919 -2777.3359 1861.2877 9.923822 17.188564 10.911763 -4214.7779 6995.1472 5169.0064 -2188.489 6157.0955 533.65478
7500 528.5729 3161.4629 -2779.0851 1855.7946 9.9140499 17.171639 10.901018 2935.365 -2873.1363 9422.1601 771.1885 -4360.9131 4939.8209
8000 533.77283 4534.849 -2777.6538 1858.4772 9.9188246 17.179909 10.906268 -1187.9433 15739.396 -946.90551 -5187.8588 2446.5059 8079.2032
8500 518.71765 1108.9877 -2777.7019 1866.6125 9.9332765 17.20494 10.922159 8720.4976 -8234.9325 2841.3979 5148.5004 -2125.3524 -4127.7468
9000 536.71495 -496.88283 -2778.0262 1877.7099 9.9529227 17.238968 10.943761 -3481.5874 -4611.6246 6602.5634 -2788.5111 -13323.148 4338.813
9500 527.06773 -236.09043 -2778.1125 1895.9227 9.9849986 17.294525 10.97903 -12233.409 7578.0514 3947.0863 -6399.0254 995.22838 8590.7109
10000 526.77335 -4480.6866 -2777.7171 1886.8998 9.9691335 17.267046 10.961585 -3139.961 1336.993 -11639.092 13496.371 -11543.676 -6180.9262
Loop time of 8.86837 on 1 procs for 10000 steps with 144 atoms
Performance: 97.425 ns/day, 0.246 hours/ns, 1127.603 timesteps/s
99.8% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 6.1503 | 6.1503 | 6.1503 | 0.0 | 69.35
Kspace | 1.1522 | 1.1522 | 1.1522 | 0.0 | 12.99
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.11705 | 0.11705 | 0.11705 | 0.0 | 1.32
Output | 0.00035834 | 0.00035834 | 0.00035834 | 0.0 | 0.00
Modify | 1.4245 | 1.4245 | 1.4245 | 0.0 | 16.06
Other | | 0.02397 | | | 0.27
Nlocal: 144 ave 144 max 144 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 3804 ave 3804 max 3804 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 41952 ave 41952 max 41952 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 41952
Ave neighs/atom = 291.333
Neighbor list builds = 0
Dangerous builds = 0
unfix quartz_qtb
unfix scapegoat_qtb
Total wall time: 0:00:08

View File

@ -0,0 +1,152 @@
LAMMPS (15 Jun 2020)
using 1 OpenMP thread(s) per MPI task
## This script first constructs an alpha quartz structure of a given size. It then uses fix qtb to equilibrate the computational cell to the specified temperature and pressure.
variable x_rep equal 2 #x-direction replication number
variable y_rep equal 2 #y-direction replication number
variable z_rep equal 2 #z-direction replication number
variable cut_off equal 10.0 #Cut-off distance for the Buckingham term (Angstrom in metal units)
variable pressure equal 1.03125 #Initial state pressure (bar in metal units)
variable temperature equal 300.0 #Initial state quantum temperature (K in metal units)
variable delta_t equal 1.0e-3 #MD timestep length (ps in metal units)
variable damp_qtb equal 1 #1/gamma where gamma is the friction coefficient in quantum thermal bath (ps in metal units)
## This part defines units, alpha-quartz crystal, and atomic information
#General
units metal
dimension 3
boundary p p p
atom_style charge
#Lattice
lattice custom 1.0 a1 4.916000 0.000000 0.000000 a2 -2.45800 4.257381 0.000000 a3 0.000000 0.000000 5.405400 basis 0.469700 0.000000 0.000000 basis 0.000000 0.469700 0.666667 basis 0.530300 0.530300 0.333333 basis 0.413500 0.266900 0.119100 basis 0.266900 0.413500 0.547567 basis 0.733100 0.146600 0.785767 basis 0.586500 0.853400 0.214233 basis 0.853400 0.586500 0.452433 basis 0.146600 0.733100 0.880900 #American Mineralogist 65 920 1980 (Space Group 154)
Lattice spacing in x,y,z = 7.374 4.25738 5.4054
#Computational Cell
region orthorhombic_unit_cell block 0 4.916000 0 8.514762 0 5.405400 units box
create_box 2 orthorhombic_unit_cell
Created orthogonal box = (0.0 0.0 0.0) to (4.916 8.514762 5.4054)
1 by 2 by 2 MPI processor grid
create_atoms 1 box basis 1 1 basis 2 1 basis 3 1 basis 4 2 basis 5 2 basis 6 2 basis 7 2 basis 8 2 basis 9 2
Created 18 atoms
create_atoms CPU = 0.000 seconds
replicate ${x_rep} ${y_rep} ${z_rep}
replicate 2 ${y_rep} ${z_rep}
replicate 2 2 ${z_rep}
replicate 2 2 2
orthogonal box = (0.0 0.0 0.0) to (9.832 17.029524 10.8108)
1 by 2 by 2 MPI processor grid
144 atoms
replicate CPU = 0.000231981 secs
#Atomic Information
mass 1 28.085500
mass 2 15.999400
set type 1 charge +2.4
48 settings made for charge
set type 2 charge -1.2
96 settings made for charge
## This part implements the BKS pair potential with a cut-off distance for the Buckingham term. Long range Coulomb interactions are evaluated with the pppm method.
#Pair Potentials
pair_style buck/coul/long ${cut_off} #BKS interaction, PRL 64 1955 (1990)
pair_style buck/coul/long 10
pair_coeff 1 1 0.0 1.0 0.0
pair_coeff 1 2 18003.757200 0.205205 133.538100
pair_coeff 2 2 1388.773000 0.362319 175.000000
pair_modify shift yes
kspace_style pppm 1.0e-4
#Neighbor style
neighbor 2.0 bin
neigh_modify check yes every 1 delay 0 page 100000 one 2000
## This part equilibrates your crystal to a pressure of ${pressure}(unit pressure) and a temperature of ${temperature}(unit temperatureture) with quantum nuclear effects
variable p_damp equal ${delta_t}*1000 #Recommended pressure damping parameter in fix nph
variable p_damp equal 0.001*1000
fix scapegoat_qtb all nph iso ${pressure} ${pressure} ${p_damp} ptemp ${temperature} #NPH does the time integration
fix scapegoat_qtb all nph iso 1.03125 ${pressure} ${p_damp} ptemp ${temperature}
fix scapegoat_qtb all nph iso 1.03125 1.03125 ${p_damp} ptemp ${temperature}
fix scapegoat_qtb all nph iso 1.03125 1.03125 1 ptemp ${temperature}
fix scapegoat_qtb all nph iso 1.03125 1.03125 1 ptemp 300
fix quartz_qtb all qtb temp ${temperature} damp ${damp_qtb} seed 35082 f_max 120.00 N_f 100 #Change f_max (THz) if your Debye frequency is higher
fix quartz_qtb all qtb temp 300 damp ${damp_qtb} seed 35082 f_max 120.00 N_f 100
fix quartz_qtb all qtb temp 300 damp 1 seed 35082 f_max 120.00 N_f 100
thermo_style custom step temp press etotal vol lx ly lz pxx pyy pzz pxy pyz pxz
thermo 500
run 10000 # 20 ps
PPPM initialization ...
using 12-bit tables for long-range coulomb (src/kspace.cpp:332)
G vector (1/distance) = 0.307414
grid = 9 15 10
stencil order = 5
estimated absolute RMS force accuracy = 0.000822922
estimated relative force accuracy = 5.71487e-05
using double precision FFTW3
3d grid and FFT values/proc = 2688 405
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 12
ghost atom cutoff = 12
binsize = 6, bins = 2 3 2
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair buck/coul/long, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
Per MPI rank memory allocation (min/avg/max) = 79.08 | 79.08 | 79.09 Mbytes
Step Temp Press TotEng Volume Lx Ly Lz Pxx Pyy Pzz Pxy Pyz Pxz
0 0 -34025.794 -2793.6041 1810.0985 9.832 17.029524 10.8108 -37478.502 -37477.413 -27121.466 -5.9958847e-10 1.3388978 7.2750373e-10
500 222.04947 -50221.579 -2787.6677 1851.5661 9.9065143 17.158586 10.892732 -61493.697 -53512.432 -35658.61 4973.9976 7095.5481 -2041.6341
1000 361.92367 -14345.85 -2783.1509 1861.0579 9.9234137 17.187857 10.911314 -4145.4149 -28701.195 -10190.939 7896.3934 -3901.2874 -490.57107
1500 457.97039 -4934.8727 -2779.8321 1860.2254 9.9219337 17.185294 10.909687 -3680.7192 -3045.0707 -8078.8283 456.70383 -4705.3346 -2175.8144
2000 523.52684 -9516.755 -2778.5181 1866.3577 9.9328244 17.204157 10.921662 -11042.489 -7777.5634 -9730.2127 -2016.3336 6027.001 -4150.3656
2500 489.58881 -4968.5157 -2777.3948 1864.0745 9.9287723 17.197139 10.917206 -13652.344 -2823.514 1570.3111 -7481.9537 -1150.3548 10502.368
3000 559.52782 -2882.7076 -2777.5527 1883.2223 9.9626528 17.255821 10.95446 3061.1755 -10570.656 -1138.642 -12045.354 -856.20951 16292.443
3500 521.67929 9974.5929 -2776.7752 1880.5936 9.9580152 17.247789 10.94936 15358.559 7855.8683 6709.3509 7292.9372 -9848.9204 -523.61056
4000 497.90872 -2012.9259 -2776.5554 1862.1703 9.9253904 17.191281 10.913488 -1154.5123 4270.0275 -9154.2927 971.94826 -10157.618 4694.0509
4500 533.64016 -7218.9278 -2775.8789 1883.3041 9.962797 17.256071 10.954618 -18299.547 -5497.566 2140.3296 -1335.6063 -10353.21 5703.7506
5000 551.61416 1590.9702 -2777.6093 1866.2047 9.9325531 17.203687 10.921363 -4600.02 6535.3 2837.6306 3412.3383 9492.18 1017.5742
5500 499.36075 188.82067 -2777.9872 1863.2925 9.9273838 17.194734 10.91568 -3238.914 1143.013 2662.363 4193.7623 -11565.423 2575.9361
6000 478.563 4064.8319 -2778.946 1867.7185 9.935238 17.208337 10.924316 1947.7246 3346.7411 6900.0301 -6339.9554 4133.6942 -4555.406
6500 512.63865 10227.461 -2778.8476 1855.5323 9.9135828 17.17083 10.900505 7423.8967 7558.2024 15700.285 -621.4585 -2620.4837 -3256.7524
7000 489.9889 13037.303 -2778.8793 1856.2469 9.9148553 17.173034 10.901904 10690.345 16770.786 11650.779 -4056.5527 -5023.8847 469.21909
7500 495.52187 1320.5068 -2778.1189 1871.7467 9.9423755 17.2207 10.932164 1978.2905 738.78041 1244.4496 1826.0923 -7829.3563 1873.2713
8000 474.60945 -4203.2068 -2778.8915 1866.5966 9.9332482 17.204891 10.922128 -1480.6896 -12516.261 1387.3306 2731.4462 -1292.9741 10743.939
8500 473.16225 -6266.1992 -2778.594 1872.9075 9.9444304 17.224259 10.934423 -12680.492 -2832.6603 -3285.4455 7226.9632 3762.6841 -5834.9064
9000 486.6579 2843.7947 -2778.0388 1877.3735 9.9523282 17.237939 10.943107 805.23659 6213.7247 1512.4228 2685.2063 -3517.5266 -17054.035
9500 549.35112 -1028.3899 -2776.8124 1880.7965 9.9583733 17.248409 10.949754 -1817.8413 2754.8459 -4022.1743 -3101.1463 8397.2345 -8608.1342
10000 562.27081 12885.53 -2775.7435 1850.2864 9.9042316 17.154633 10.890222 15758.218 9989.5121 12908.859 -25.724137 -16691.374 267.85371
Loop time of 3.80648 on 4 procs for 10000 steps with 144 atoms
Performance: 226.981 ns/day, 0.106 hours/ns, 2627.100 timesteps/s
94.8% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 1.3526 | 1.6581 | 1.9634 | 21.3 | 43.56
Kspace | 0.92143 | 1.2222 | 1.5232 | 24.7 | 32.11
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.31886 | 0.32256 | 0.32604 | 0.5 | 8.47
Output | 0.0003643 | 0.00083923 | 0.0022533 | 0.0 | 0.02
Modify | 0.39166 | 0.45985 | 0.52607 | 8.9 | 12.08
Other | | 0.143 | | | 3.76
Nlocal: 36 ave 36 max 36 min
Histogram: 4 0 0 0 0 0 0 0 0 0
Nghost: 2614 ave 2614 max 2614 min
Histogram: 4 0 0 0 0 0 0 0 0 0
Neighs: 10488 ave 12570 max 8406 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Total # of neighbors = 41952
Ave neighs/atom = 291.333
Neighbor list builds = 0
Dangerous builds = 0
unfix quartz_qtb
unfix scapegoat_qtb
Total wall time: 0:00:03

View File

@ -0,0 +1,183 @@
Reactive MD-force field
39 ! Number of general parameters
50.0000 !Overcoordination parameter
9.5469 !Overcoordination parameter
127.8302 !Valency angle conjugation parameter
3.0000 !Triple bond stabilisation parameter
6.5000 !Triple bond stabilisation parameter
0.0000 !C2-correction
1.0496 !Undercoordination parameter
9.0000 !Triple bond stabilisation parameter
11.5054 !Undercoordination parameter
13.4059 !Undercoordination parameter
0.0000 !Triple bond stabilization energy
0.0000 !Lower Taper-radius
10.0000 !Upper Taper-radius
2.8793 !Not used
33.8667 !Valency undercoordination
7.0994 !Valency angle/lone pair parameter
1.0563 !Valency angle
2.0384 !Valency angle parameter
6.1431 !Not used
6.9290 !Double bond/angle parameter
0.3989 !Double bond/angle parameter: overcoord
3.9954 !Double bond/angle parameter: overcoord
-2.4837 !Not used
5.7796 !Torsion/BO parameter
10.0000 !Torsion overcoordination
1.9487 !Torsion overcoordination
-1.2327 !Conjugation 0 (not used)
2.1645 !Conjugation
1.5591 !vdWaals shielding
0.1000 !Cutoff for bond order (*100)
2.0038 !Valency angle conjugation parameter
0.6121 !Overcoordination parameter
1.2172 !Overcoordination parameter
1.8512 !Valency/lone pair parameter
0.5000 !Not used
20.0000 !Not used
5.0000 !Molecular energy (not used)
0.0000 !Molecular energy (not used)
3.6942 !Valency angle conjugation parameter
5 ! Nr of atoms; cov.r; valency;a.m;Rvdw;Evdw;gammaEEM;cov.r2;#
alfa;gammavdW;valency;Eunder;Eover;chiEEM;etaEEM;n.u.
cov r3;Elp;Heat inc.;n.u.;n.u.;n.u.;n.u.
ov/un;val1;n.u.;val3,vval4
C 1.3763 4.0000 12.0000 1.8857 0.1818 0.8712 1.2596 4.0000
9.5928 2.0784 4.0000 22.6732 79.5548 5.7254 6.9235 0.0000
1.2065 0.0000 -0.8579 4.9417 28.3475 11.9957 0.8563 0.0000
-2.8846 4.1590 1.0564 4.0000 2.9663 0.0000 0.0000 0.0000
H 0.6646 1.0000 1.0080 1.6030 0.0600 0.7625 -0.1000 1.0000
9.3951 4.4187 1.0000 0.0000 121.1250 3.8196 9.8832 1.0000
-0.1000 0.0000 -0.1339 3.5803 2.8733 1.0000 1.0698 0.0000
-13.0615 3.0626 1.0338 1.0000 2.8793 0.0000 0.0000 0.0000
O 1.2699 2.0000 15.9990 1.9741 0.0880 1.0804 1.0624 6.0000
10.2186 7.7719 4.0000 27.3264 116.0768 8.5000 7.8386 2.0000
0.9446 8.6170 -1.2371 17.0845 3.7082 0.5350 0.9745 0.0000
-3.1456 2.6656 1.0493 4.0000 2.9225 0.0000 0.0000 0.0000
N 1.2226 3.0000 14.0000 1.9324 0.1376 0.8596 1.1839 5.0000
10.0667 7.8431 4.0000 32.5000 100.0000 6.8418 6.3404 2.0000
1.0497 14.5853 -1.1222 2.0637 3.2584 3.1136 0.9745 0.0000
-4.2059 2.6491 1.0183 4.0000 2.8793 0.0000 0.0000 0.0000
S 1.9405 2.0000 32.0600 2.0677 0.2099 1.0336 1.5479 6.0000
9.9575 4.9055 4.0000 52.9998 112.1416 6.5000 8.2545 2.0000
1.4601 9.7177 -2.3700 5.7487 23.2859 12.7147 0.9745 0.0000
-11.0000 2.7466 1.0338 4.0000 2.8793 0.0000 0.0000 0.0000
15 ! Nr of bonds; Edis1;LPpen;n.u.;pbe1;pbo5;13corr;pbo6
pbe2;pbo3;pbo4;n.u.;pbo1;pbo2;ovcorr
1 1 145.4070 103.0681 73.7841 0.2176 -0.7816 1.0000 28.4167 0.3217
0.1111 -0.1940 8.6733 1.0000 -0.0994 5.9724 1.0000 0.0000
1 2 167.1752 0.0000 0.0000 -0.4421 0.0000 1.0000 6.0000 0.5969
17.4194 1.0000 0.0000 1.0000 -0.0099 8.5445 0.0000 0.0000
2 2 188.1606 0.0000 0.0000 -0.3140 0.0000 1.0000 6.0000 0.6816
8.6247 1.0000 0.0000 1.0000 -0.0183 5.7082 0.0000 0.0000
1 3 171.0470 67.2480 130.3792 0.3600 -0.1696 1.0000 12.0338 0.3796
0.3647 -0.2660 7.4396 1.0000 -0.1661 5.0637 0.0000 0.0000
3 3 90.2465 160.9645 40.0000 0.9950 -0.2435 1.0000 28.1614 0.9704
0.8145 -0.1850 7.5281 1.0000 -0.1283 6.2396 1.0000 0.0000
1 4 134.9992 139.6314 78.5681 0.0420 -0.1370 1.0000 23.6247 0.2415
0.1522 -0.3161 7.0000 1.0000 -0.1301 5.4980 1.0000 0.0000
3 4 127.7074 177.1058 40.0000 0.4561 -0.1481 1.0000 31.4801 0.2000
0.8968 -0.3555 7.0000 1.0000 -0.1219 7.0000 1.0000 0.0000
4 4 151.9142 87.1928 151.4761 0.4280 -0.1001 1.0000 12.3631 0.6229
0.1721 -0.1614 12.1345 1.0000 -0.0882 5.3056 1.0000 0.0000
2 3 216.6018 0.0000 0.0000 -0.4201 0.0000 1.0000 6.0000 0.9143
4.7737 1.0000 0.0000 1.0000 -0.0591 5.9451 0.0000 0.0000
2 4 223.1853 0.0000 0.0000 -0.4661 0.0000 1.0000 6.0000 0.5178
7.8731 1.0000 0.0000 1.0000 -0.0306 6.1506 0.0000 0.0000
1 5 128.9942 74.5848 55.2528 0.1035 -0.5211 1.0000 18.9617 0.6000
0.2949 -0.2398 8.1175 1.0000 -0.1029 5.6731 1.0000 0.0000
2 5 151.5159 0.0000 0.0000 -0.4721 0.0000 1.0000 6.0000 0.6000
9.4366 1.0000 0.0000 1.0000 -0.0290 7.0050 1.0000 0.0000
3 5 0.0000 0.0000 0.0000 0.5563 -0.4038 1.0000 49.5611 0.6000
0.4259 -0.4577 12.7569 1.0000 -0.1100 7.1145 1.0000 0.0000
4 5 0.0000 0.0000 0.0000 0.4438 -0.2034 1.0000 40.3399 0.6000
0.3296 -0.3153 9.1227 1.0000 -0.1805 5.6864 1.0000 0.0000
5 5 96.1871 93.7006 68.6860 0.0955 -0.4781 1.0000 17.8574 0.6000
0.2723 -0.2373 9.7875 1.0000 -0.0950 6.4757 1.0000 0.0000
6 ! Nr of off-diagonal terms; Ediss;Ro;gamma;rsigma;rpi;rpi2
1 2 0.0455 1.7218 10.4236 1.0379 -1.0000 -1.0000
2 3 0.0469 1.9185 10.3707 0.9406 -1.0000 -1.0000
2 4 0.0999 1.8372 9.6539 0.9692 -1.0000 -1.0000
1 3 0.1186 1.9820 9.5927 1.2936 1.1203 1.0805
1 4 0.1486 1.8922 9.7989 1.3746 1.2091 1.1427
3 4 0.1051 2.0060 10.0691 1.3307 1.1034 1.0060
50 ! Nr of angles;at1;at2;at3;Thetao,o;ka;kb;pv1;pv2
1 1 1 70.0265 13.6338 2.1884 0.0000 0.1676 26.3587 1.0400
1 1 2 69.7786 10.3544 8.4326 0.0000 0.1153 0.0000 1.0400
2 1 2 74.6020 11.8629 2.9294 0.0000 0.1367 0.0000 1.0400
1 2 2 0.0000 0.0000 6.0000 0.0000 0.0000 0.0000 1.0400
1 2 1 0.0000 3.4110 7.7350 0.0000 0.0000 0.0000 1.0400
2 2 2 0.0000 27.9213 5.8635 0.0000 0.0000 0.0000 1.0400
1 1 3 72.9588 16.7105 3.5244 0.0000 1.1127 0.0000 1.1880
3 1 3 80.0708 45.0000 2.1487 0.0000 1.1127 -35.0000 1.1880
1 1 4 61.5055 45.0000 1.2242 0.0000 1.1127 0.0000 1.1880
3 1 4 71.9345 45.0000 1.5052 0.0000 1.1127 0.0000 1.1880
4 1 4 51.3604 45.0000 0.6846 0.0000 1.1127 0.0000 1.1880
2 1 3 66.6150 13.6403 3.8212 0.0000 0.0755 0.0000 1.0500
2 1 4 68.9632 16.3575 3.1449 0.0000 0.0755 0.0000 1.0500
1 2 4 0.0000 0.0019 6.3000 0.0000 0.0000 0.0000 1.0400
1 3 1 79.1091 45.0000 0.7067 0.0000 0.6142 0.0000 1.0783
1 3 3 83.7151 42.6867 0.9699 0.0000 0.6142 0.0000 1.0783
1 3 4 79.5876 45.0000 1.1761 0.0000 0.6142 0.0000 1.0783
3 3 3 80.0108 38.3716 1.1572 -38.4200 0.6142 0.0000 1.0783
3 3 4 81.5614 19.8012 3.9968 0.0000 0.6142 0.0000 1.0783
4 3 4 85.3564 36.5858 1.7504 0.0000 0.6142 0.0000 1.0783
1 3 2 78.1533 44.7226 1.3136 0.0000 0.1218 0.0000 1.0500
2 3 3 84.1057 9.6413 7.5000 0.0000 0.1218 0.0000 1.0500
2 3 4 79.4629 44.0409 2.2959 0.0000 0.1218 0.0000 1.0500
2 3 2 79.2954 26.3838 2.2044 0.0000 0.1218 0.0000 1.0500
1 4 1 66.1477 22.9891 1.5923 0.0000 1.6777 0.0000 1.0500
1 4 3 91.9273 38.0207 0.5387 0.0000 1.6777 0.0000 1.0500
1 4 4 92.6933 9.9708 1.6094 0.0000 1.6777 0.0000 1.0500
3 4 3 73.4749 42.7640 1.7325 -17.5007 1.6777 0.0000 1.0500
3 4 4 73.9183 44.8857 1.1980 -0.9193 1.6777 0.0000 1.0500
4 4 4 74.0572 15.4709 5.4220 0.0000 1.6777 0.0000 1.0500
1 4 2 72.7016 33.4153 1.0224 0.0000 0.0222 0.0000 1.0500
2 4 3 82.4368 44.1900 1.9273 0.0000 0.0222 0.0000 1.0500
2 4 4 82.6883 39.9831 1.1916 0.0000 0.0222 0.0000 1.0500
2 4 2 71.2183 14.4528 3.6870 0.0000 0.0222 0.0000 1.0500
1 2 3 0.0000 0.0019 6.0000 0.0000 0.0000 0.0000 1.0400
1 2 4 0.0000 0.0019 6.0000 0.0000 0.0000 0.0000 1.0400
1 2 5 0.0000 0.0019 6.0000 0.0000 0.0000 0.0000 1.0400
3 2 3 0.0000 0.0019 6.0000 0.0000 0.0000 0.0000 1.0400
3 2 4 0.0000 0.0019 6.0000 0.0000 0.0000 0.0000 1.0400
4 2 4 0.0000 0.0019 6.0000 0.0000 0.0000 0.0000 1.0400
2 2 3 0.0000 0.0019 6.0000 0.0000 0.0000 0.0000 1.0400
2 2 4 0.0000 0.0019 6.0000 0.0000 0.0000 0.0000 1.0400
1 1 5 74.9397 25.0560 1.8787 0.1463 0.0559 0.0000 1.0400
1 5 1 86.9521 36.9951 2.0903 0.1463 0.0559 0.0000 1.0400
2 1 5 74.9397 25.0560 1.8787 0.0000 0.0000 0.0000 1.0400
1 5 2 86.1791 36.9951 2.0903 0.0000 0.0000 0.0000 1.0400
1 5 5 85.3644 36.9951 2.0903 0.1463 0.0559 0.0000 1.0400
2 5 2 93.1959 36.9951 2.0903 0.0000 0.0000 0.0000 1.0400
2 5 5 84.3331 36.9951 2.0903 0.0000 0.0000 0.0000 1.0400
2 2 5 0.0000 0.0019 6.0000 0.0000 0.0000 0.0000 1.0400
17 ! Nr of torsions;at1;at2;at3;at4;;V1;V2;V3;V2(BO);vconj;n.u;n
1 1 1 1 0.0000 23.2168 0.1811 -4.6220 -1.9387 0.0000 0.0000
1 1 1 2 0.0000 45.7984 0.3590 -5.7106 -2.9459 0.0000 0.0000
2 1 1 2 0.0000 44.6445 0.3486 -5.1725 -0.8717 0.0000 0.0000
0 1 2 0 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
0 2 2 0 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
0 1 3 0 5.0520 16.7344 0.5590 -3.0181 -2.0000 0.0000 0.0000
0 2 3 0 0.0000 0.1000 0.0200 -2.5415 0.0000 0.0000 0.0000
0 3 3 0 0.0115 68.9706 0.8253 -28.4693 0.0000 0.0000 0.0000
0 1 4 0 -4.0616 66.2036 0.3855 -4.4414 -2.0000 0.0000 0.0000
0 2 4 0 0.0000 0.1000 0.0200 -2.5415 0.0000 0.0000 0.0000
0 3 4 0 1.1130 14.8049 0.0231 -10.7175 -2.0000 0.0000 0.0000
0 4 4 0 -0.0851 37.4200 0.0107 -3.5209 -2.0000 0.0000 0.0000
0 1 1 0 0.0000 0.9305 0.0000 -24.2568 0.0000 0.0000 0.0000
4 1 4 4 -3.6064 43.6430 0.0004 -11.5507 -2.0000 0.0000 0.0000
0 1 5 0 3.3423 30.3435 0.0365 -2.7171 0.0000 0.0000 0.0000
0 5 5 0 -0.0555 -42.7738 0.1515 -2.2056 0.0000 0.0000 0.0000
0 2 5 0 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
9 ! Nr of hydrogen bonds;at1;at2;at3;Rhb;Dehb;vhb1
3 2 3 2.0431 -6.6813 3.5000 1.7295
3 2 4 1.6740 -10.9581 3.5000 1.7295
4 2 3 1.4889 -9.6465 3.5000 1.7295
4 2 4 1.8324 -8.0074 3.5000 1.7295
3 2 5 2.6644 -3.9547 3.5000 1.7295
4 2 5 4.0476 -5.7038 3.5000 1.7295
5 2 3 2.1126 -4.5790 3.5000 1.7295
5 2 4 2.2066 -5.7038 3.5000 1.7295
5 2 5 1.9461 -4.0000 3.5000 1.7295

View File

@ -0,0 +1,33 @@
## This script first uses fix qtb to equilibrate liquid methane to an initial state with quantum nuclear correction and then simulate shock induced chemical reactions through the quantum thermal bath multi-scale shock technique
#The default system size may take a while to run you can change to a smaller size
variable x_rep equal 5 #x-direction replication number
variable y_rep equal 5 #y-direction replication number
variable z_rep equal 10 #z-direction replication number
variable temperature equal 110.0 #Target quantum temperature (K in real units)
variable delta_t equal 0.25 #MD timestep length (fs in real units)
variable damp_qtb equal 200 #1/gamma where gamma is the friction coefficient in quantum thermal bath (fs in real units)
variable v_msst equal 0.122 #Shock velocity (Angstrom/fs in metal units)
variable q_msst equal 25.0 #Box mass-like parameter in the MSST (mass^2/length^4, where mass=grams/mole and length=Angstrom in real units)
variable mu_msst equal 0.9 #Artificial viscosity in the MSST (mass/length/time, where mass=grams/mole, length=Angstrom and time=fs in real units)
variable tscale_msst equal 0.01 #Temperature reduction parameter in the MSST (unitless)
variable eta_qbmsst equal 1.0 #Coupling constant between the shock and the quantum thermal bath (unitless constant)
##The included part first constructs a liquid methane structure of a given size. It then uses fix qtb to equilibrate the computational cell to the specified temperature and pressure.
include methane_qtb.mod
##Shock compression with quantum nuclear corrections
reset_timestep 0
fix shock all qbmsst z ${v_msst} q ${q_msst} mu ${mu_msst} tscale ${tscale_msst} damp ${damp_qtb} f_max 0.3 N_f 50 seed 35082 eta ${eta_qbmsst} beta 400 T_init ${temperature}
fix_modify shock energy yes
variable dhug equal f_shock[1]
variable dray equal f_shock[2]
variable lgr_vel equal f_shock[3]
variable lgr_pos equal f_shock[4]
variable T_qm equal f_shock[5] #Temperature with quantum nuclear correction
thermo_style custom step v_T_qm press etotal vol lx ly lz pzz v_dhug v_dray v_lgr_vel v_lgr_pos
thermo 100
timestep ${delta_t}
restart 1000 restart
run 5000

View File

@ -0,0 +1,64 @@
## This script first constructs a liquid methane structure of a given size. It then uses fix qtb to equilibrate the computational cell to the specified temperature and pressure.
## This part defines units, methane structure, and atomic information
#General
units real
dimension 3
boundary p p p
atom_style charge
#Lattice
lattice custom 1.0 &
a1 3.9783624 0 0 &
a2 0 3.9783624 0 &
a3 0 0 3.9783624 &
&
basis 0.5 0.5 0.5 &
basis 0.663 0.663 0.663 &
basis 0.337 0.337 0.663 &
basis 0.663 0.337 0.337 &
basis 0.337 0.663 0.337
#Computational Cell
region simbox block 0 3.9783624 0 3.9783624 0 3.9783624 units box
create_box 2 simbox
create_atoms 1 box &
basis 1 1 &
basis 2 2 &
basis 3 2 &
basis 4 2 &
basis 5 2
replicate ${x_rep} ${y_rep} ${z_rep}
#Atomic Information
mass 1 12.011150
mass 2 1.007970
## This part defines the reax pair potential in methane, force field coefficients are specified in "ffield.reax"
#Pair Potentials
pair_style reax/c NULL
pair_coeff * * ffield.reax C H
fix 0 all qeq/reax 1 0.0 10.0 1.0e-6 reax/c
#Neighbor Style
neighbor 2.5 bin
neigh_modify every 10 delay 0 check no
## This part equilibrates liquid methane to a temperature of ${temperature}(unit temperatureture) with quantum nuclear effects
#Initialization
velocity all create ${temperature} 93 dist gaussian sum no mom yes rot yes loop all
#Setup output
thermo_style custom step temp press etotal vol
thermo 100
#Colored thermal bath
fix scapegoat_qtb all nve #NVE does the time integration
fix methane_qtb all qtb temp ${temperature} damp ${damp_qtb} seed 35082 f_max 0.3 N_f 50 #Change f_max if your Debye frequency is higher
timestep ${delta_t}
run 2000 #500 fs
unfix methane_qtb
unfix scapegoat_qtb

View File

@ -0,0 +1,183 @@
Reactive MD-force field
39 ! Number of general parameters
50.0000 !Overcoordination parameter
9.5469 !Overcoordination parameter
127.8302 !Valency angle conjugation parameter
3.0000 !Triple bond stabilisation parameter
6.5000 !Triple bond stabilisation parameter
0.0000 !C2-correction
1.0496 !Undercoordination parameter
9.0000 !Triple bond stabilisation parameter
11.5054 !Undercoordination parameter
13.4059 !Undercoordination parameter
0.0000 !Triple bond stabilization energy
0.0000 !Lower Taper-radius
10.0000 !Upper Taper-radius
2.8793 !Not used
33.8667 !Valency undercoordination
7.0994 !Valency angle/lone pair parameter
1.0563 !Valency angle
2.0384 !Valency angle parameter
6.1431 !Not used
6.9290 !Double bond/angle parameter
0.3989 !Double bond/angle parameter: overcoord
3.9954 !Double bond/angle parameter: overcoord
-2.4837 !Not used
5.7796 !Torsion/BO parameter
10.0000 !Torsion overcoordination
1.9487 !Torsion overcoordination
-1.2327 !Conjugation 0 (not used)
2.1645 !Conjugation
1.5591 !vdWaals shielding
0.1000 !Cutoff for bond order (*100)
2.0038 !Valency angle conjugation parameter
0.6121 !Overcoordination parameter
1.2172 !Overcoordination parameter
1.8512 !Valency/lone pair parameter
0.5000 !Not used
20.0000 !Not used
5.0000 !Molecular energy (not used)
0.0000 !Molecular energy (not used)
3.6942 !Valency angle conjugation parameter
5 ! Nr of atoms; cov.r; valency;a.m;Rvdw;Evdw;gammaEEM;cov.r2;#
alfa;gammavdW;valency;Eunder;Eover;chiEEM;etaEEM;n.u.
cov r3;Elp;Heat inc.;n.u.;n.u.;n.u.;n.u.
ov/un;val1;n.u.;val3,vval4
C 1.3763 4.0000 12.0000 1.8857 0.1818 0.8712 1.2596 4.0000
9.5928 2.0784 4.0000 22.6732 79.5548 5.7254 6.9235 0.0000
1.2065 0.0000 -0.8579 4.9417 28.3475 11.9957 0.8563 0.0000
-2.8846 4.1590 1.0564 4.0000 2.9663 0.0000 0.0000 0.0000
H 0.6646 1.0000 1.0080 1.6030 0.0600 0.7625 -0.1000 1.0000
9.3951 4.4187 1.0000 0.0000 121.1250 3.8196 9.8832 1.0000
-0.1000 0.0000 -0.1339 3.5803 2.8733 1.0000 1.0698 0.0000
-13.0615 3.0626 1.0338 1.0000 2.8793 0.0000 0.0000 0.0000
O 1.2699 2.0000 15.9990 1.9741 0.0880 1.0804 1.0624 6.0000
10.2186 7.7719 4.0000 27.3264 116.0768 8.5000 7.8386 2.0000
0.9446 8.6170 -1.2371 17.0845 3.7082 0.5350 0.9745 0.0000
-3.1456 2.6656 1.0493 4.0000 2.9225 0.0000 0.0000 0.0000
N 1.2226 3.0000 14.0000 1.9324 0.1376 0.8596 1.1839 5.0000
10.0667 7.8431 4.0000 32.5000 100.0000 6.8418 6.3404 2.0000
1.0497 14.5853 -1.1222 2.0637 3.2584 3.1136 0.9745 0.0000
-4.2059 2.6491 1.0183 4.0000 2.8793 0.0000 0.0000 0.0000
S 1.9405 2.0000 32.0600 2.0677 0.2099 1.0336 1.5479 6.0000
9.9575 4.9055 4.0000 52.9998 112.1416 6.5000 8.2545 2.0000
1.4601 9.7177 -2.3700 5.7487 23.2859 12.7147 0.9745 0.0000
-11.0000 2.7466 1.0338 4.0000 2.8793 0.0000 0.0000 0.0000
15 ! Nr of bonds; Edis1;LPpen;n.u.;pbe1;pbo5;13corr;pbo6
pbe2;pbo3;pbo4;n.u.;pbo1;pbo2;ovcorr
1 1 145.4070 103.0681 73.7841 0.2176 -0.7816 1.0000 28.4167 0.3217
0.1111 -0.1940 8.6733 1.0000 -0.0994 5.9724 1.0000 0.0000
1 2 167.1752 0.0000 0.0000 -0.4421 0.0000 1.0000 6.0000 0.5969
17.4194 1.0000 0.0000 1.0000 -0.0099 8.5445 0.0000 0.0000
2 2 188.1606 0.0000 0.0000 -0.3140 0.0000 1.0000 6.0000 0.6816
8.6247 1.0000 0.0000 1.0000 -0.0183 5.7082 0.0000 0.0000
1 3 171.0470 67.2480 130.3792 0.3600 -0.1696 1.0000 12.0338 0.3796
0.3647 -0.2660 7.4396 1.0000 -0.1661 5.0637 0.0000 0.0000
3 3 90.2465 160.9645 40.0000 0.9950 -0.2435 1.0000 28.1614 0.9704
0.8145 -0.1850 7.5281 1.0000 -0.1283 6.2396 1.0000 0.0000
1 4 134.9992 139.6314 78.5681 0.0420 -0.1370 1.0000 23.6247 0.2415
0.1522 -0.3161 7.0000 1.0000 -0.1301 5.4980 1.0000 0.0000
3 4 127.7074 177.1058 40.0000 0.4561 -0.1481 1.0000 31.4801 0.2000
0.8968 -0.3555 7.0000 1.0000 -0.1219 7.0000 1.0000 0.0000
4 4 151.9142 87.1928 151.4761 0.4280 -0.1001 1.0000 12.3631 0.6229
0.1721 -0.1614 12.1345 1.0000 -0.0882 5.3056 1.0000 0.0000
2 3 216.6018 0.0000 0.0000 -0.4201 0.0000 1.0000 6.0000 0.9143
4.7737 1.0000 0.0000 1.0000 -0.0591 5.9451 0.0000 0.0000
2 4 223.1853 0.0000 0.0000 -0.4661 0.0000 1.0000 6.0000 0.5178
7.8731 1.0000 0.0000 1.0000 -0.0306 6.1506 0.0000 0.0000
1 5 128.9942 74.5848 55.2528 0.1035 -0.5211 1.0000 18.9617 0.6000
0.2949 -0.2398 8.1175 1.0000 -0.1029 5.6731 1.0000 0.0000
2 5 151.5159 0.0000 0.0000 -0.4721 0.0000 1.0000 6.0000 0.6000
9.4366 1.0000 0.0000 1.0000 -0.0290 7.0050 1.0000 0.0000
3 5 0.0000 0.0000 0.0000 0.5563 -0.4038 1.0000 49.5611 0.6000
0.4259 -0.4577 12.7569 1.0000 -0.1100 7.1145 1.0000 0.0000
4 5 0.0000 0.0000 0.0000 0.4438 -0.2034 1.0000 40.3399 0.6000
0.3296 -0.3153 9.1227 1.0000 -0.1805 5.6864 1.0000 0.0000
5 5 96.1871 93.7006 68.6860 0.0955 -0.4781 1.0000 17.8574 0.6000
0.2723 -0.2373 9.7875 1.0000 -0.0950 6.4757 1.0000 0.0000
6 ! Nr of off-diagonal terms; Ediss;Ro;gamma;rsigma;rpi;rpi2
1 2 0.0455 1.7218 10.4236 1.0379 -1.0000 -1.0000
2 3 0.0469 1.9185 10.3707 0.9406 -1.0000 -1.0000
2 4 0.0999 1.8372 9.6539 0.9692 -1.0000 -1.0000
1 3 0.1186 1.9820 9.5927 1.2936 1.1203 1.0805
1 4 0.1486 1.8922 9.7989 1.3746 1.2091 1.1427
3 4 0.1051 2.0060 10.0691 1.3307 1.1034 1.0060
50 ! Nr of angles;at1;at2;at3;Thetao,o;ka;kb;pv1;pv2
1 1 1 70.0265 13.6338 2.1884 0.0000 0.1676 26.3587 1.0400
1 1 2 69.7786 10.3544 8.4326 0.0000 0.1153 0.0000 1.0400
2 1 2 74.6020 11.8629 2.9294 0.0000 0.1367 0.0000 1.0400
1 2 2 0.0000 0.0000 6.0000 0.0000 0.0000 0.0000 1.0400
1 2 1 0.0000 3.4110 7.7350 0.0000 0.0000 0.0000 1.0400
2 2 2 0.0000 27.9213 5.8635 0.0000 0.0000 0.0000 1.0400
1 1 3 72.9588 16.7105 3.5244 0.0000 1.1127 0.0000 1.1880
3 1 3 80.0708 45.0000 2.1487 0.0000 1.1127 -35.0000 1.1880
1 1 4 61.5055 45.0000 1.2242 0.0000 1.1127 0.0000 1.1880
3 1 4 71.9345 45.0000 1.5052 0.0000 1.1127 0.0000 1.1880
4 1 4 51.3604 45.0000 0.6846 0.0000 1.1127 0.0000 1.1880
2 1 3 66.6150 13.6403 3.8212 0.0000 0.0755 0.0000 1.0500
2 1 4 68.9632 16.3575 3.1449 0.0000 0.0755 0.0000 1.0500
1 2 4 0.0000 0.0019 6.3000 0.0000 0.0000 0.0000 1.0400
1 3 1 79.1091 45.0000 0.7067 0.0000 0.6142 0.0000 1.0783
1 3 3 83.7151 42.6867 0.9699 0.0000 0.6142 0.0000 1.0783
1 3 4 79.5876 45.0000 1.1761 0.0000 0.6142 0.0000 1.0783
3 3 3 80.0108 38.3716 1.1572 -38.4200 0.6142 0.0000 1.0783
3 3 4 81.5614 19.8012 3.9968 0.0000 0.6142 0.0000 1.0783
4 3 4 85.3564 36.5858 1.7504 0.0000 0.6142 0.0000 1.0783
1 3 2 78.1533 44.7226 1.3136 0.0000 0.1218 0.0000 1.0500
2 3 3 84.1057 9.6413 7.5000 0.0000 0.1218 0.0000 1.0500
2 3 4 79.4629 44.0409 2.2959 0.0000 0.1218 0.0000 1.0500
2 3 2 79.2954 26.3838 2.2044 0.0000 0.1218 0.0000 1.0500
1 4 1 66.1477 22.9891 1.5923 0.0000 1.6777 0.0000 1.0500
1 4 3 91.9273 38.0207 0.5387 0.0000 1.6777 0.0000 1.0500
1 4 4 92.6933 9.9708 1.6094 0.0000 1.6777 0.0000 1.0500
3 4 3 73.4749 42.7640 1.7325 -17.5007 1.6777 0.0000 1.0500
3 4 4 73.9183 44.8857 1.1980 -0.9193 1.6777 0.0000 1.0500
4 4 4 74.0572 15.4709 5.4220 0.0000 1.6777 0.0000 1.0500
1 4 2 72.7016 33.4153 1.0224 0.0000 0.0222 0.0000 1.0500
2 4 3 82.4368 44.1900 1.9273 0.0000 0.0222 0.0000 1.0500
2 4 4 82.6883 39.9831 1.1916 0.0000 0.0222 0.0000 1.0500
2 4 2 71.2183 14.4528 3.6870 0.0000 0.0222 0.0000 1.0500
1 2 3 0.0000 0.0019 6.0000 0.0000 0.0000 0.0000 1.0400
1 2 4 0.0000 0.0019 6.0000 0.0000 0.0000 0.0000 1.0400
1 2 5 0.0000 0.0019 6.0000 0.0000 0.0000 0.0000 1.0400
3 2 3 0.0000 0.0019 6.0000 0.0000 0.0000 0.0000 1.0400
3 2 4 0.0000 0.0019 6.0000 0.0000 0.0000 0.0000 1.0400
4 2 4 0.0000 0.0019 6.0000 0.0000 0.0000 0.0000 1.0400
2 2 3 0.0000 0.0019 6.0000 0.0000 0.0000 0.0000 1.0400
2 2 4 0.0000 0.0019 6.0000 0.0000 0.0000 0.0000 1.0400
1 1 5 74.9397 25.0560 1.8787 0.1463 0.0559 0.0000 1.0400
1 5 1 86.9521 36.9951 2.0903 0.1463 0.0559 0.0000 1.0400
2 1 5 74.9397 25.0560 1.8787 0.0000 0.0000 0.0000 1.0400
1 5 2 86.1791 36.9951 2.0903 0.0000 0.0000 0.0000 1.0400
1 5 5 85.3644 36.9951 2.0903 0.1463 0.0559 0.0000 1.0400
2 5 2 93.1959 36.9951 2.0903 0.0000 0.0000 0.0000 1.0400
2 5 5 84.3331 36.9951 2.0903 0.0000 0.0000 0.0000 1.0400
2 2 5 0.0000 0.0019 6.0000 0.0000 0.0000 0.0000 1.0400
17 ! Nr of torsions;at1;at2;at3;at4;;V1;V2;V3;V2(BO);vconj;n.u;n
1 1 1 1 0.0000 23.2168 0.1811 -4.6220 -1.9387 0.0000 0.0000
1 1 1 2 0.0000 45.7984 0.3590 -5.7106 -2.9459 0.0000 0.0000
2 1 1 2 0.0000 44.6445 0.3486 -5.1725 -0.8717 0.0000 0.0000
0 1 2 0 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
0 2 2 0 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
0 1 3 0 5.0520 16.7344 0.5590 -3.0181 -2.0000 0.0000 0.0000
0 2 3 0 0.0000 0.1000 0.0200 -2.5415 0.0000 0.0000 0.0000
0 3 3 0 0.0115 68.9706 0.8253 -28.4693 0.0000 0.0000 0.0000
0 1 4 0 -4.0616 66.2036 0.3855 -4.4414 -2.0000 0.0000 0.0000
0 2 4 0 0.0000 0.1000 0.0200 -2.5415 0.0000 0.0000 0.0000
0 3 4 0 1.1130 14.8049 0.0231 -10.7175 -2.0000 0.0000 0.0000
0 4 4 0 -0.0851 37.4200 0.0107 -3.5209 -2.0000 0.0000 0.0000
0 1 1 0 0.0000 0.9305 0.0000 -24.2568 0.0000 0.0000 0.0000
4 1 4 4 -3.6064 43.6430 0.0004 -11.5507 -2.0000 0.0000 0.0000
0 1 5 0 3.3423 30.3435 0.0365 -2.7171 0.0000 0.0000 0.0000
0 5 5 0 -0.0555 -42.7738 0.1515 -2.2056 0.0000 0.0000 0.0000
0 2 5 0 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
9 ! Nr of hydrogen bonds;at1;at2;at3;Rhb;Dehb;vhb1
3 2 3 2.0431 -6.6813 3.5000 1.7295
3 2 4 1.6740 -10.9581 3.5000 1.7295
4 2 3 1.4889 -9.6465 3.5000 1.7295
4 2 4 1.8324 -8.0074 3.5000 1.7295
3 2 5 2.6644 -3.9547 3.5000 1.7295
4 2 5 4.0476 -5.7038 3.5000 1.7295
5 2 3 2.1126 -4.5790 3.5000 1.7295
5 2 4 2.2066 -5.7038 3.5000 1.7295
5 2 5 1.9461 -4.0000 3.5000 1.7295

View File

@ -0,0 +1,70 @@
## This script first constructs a liquid methane structure of a given size. It then uses fix qtb to equilibrate the computational cell to the specified temperature and pressure.
variable x_rep equal 2 #x-direction replication number
variable y_rep equal 2 #y-direction replication number
variable z_rep equal 2 #z-direction replication number
variable temperature equal 110.0 #Target quantum temperature (K in real units)
variable delta_t equal 0.25 #MD timestep length (fs in real units)
variable damp_qtb equal 200 #1/gamma where gamma is the friction coefficient in quantum thermal bath (fs in real units)
## This part defines units, methane structure, and atomic information
#General
units real
dimension 3
boundary p p p
atom_style charge
#Lattice
lattice custom 1.0 &
a1 3.9783624 0 0 &
a2 0 3.9783624 0 &
a3 0 0 3.9783624 &
&
basis 0.5 0.5 0.5 &
basis 0.663 0.663 0.663 &
basis 0.337 0.337 0.663 &
basis 0.663 0.337 0.337 &
basis 0.337 0.663 0.337
#Computational Cell
region simbox block 0 3.9783624 0 3.9783624 0 3.9783624 units box
create_box 2 simbox
create_atoms 1 box &
basis 1 1 &
basis 2 2 &
basis 3 2 &
basis 4 2 &
basis 5 2
replicate ${x_rep} ${y_rep} ${z_rep}
#Atomic Information
mass 1 12.011150
mass 2 1.007970
## This part defines the reax pair potential in methane, force field coefficients are specified in "ffield.reax"
#Pair Potentials
pair_style reax/c NULL
pair_coeff * * ffield.reax C H
fix 0 all qeq/reax 1 0.0 10.0 1.0e-6 reax/c
#Neighbor Style
neighbor 2.5 bin
neigh_modify every 10 delay 0 check no
## This part equilibrates liquid methane to a temperature of ${temperature}(unit temperatureture) with quantum nuclear effects
#Initialization
velocity all create ${temperature} 93 dist gaussian sum no mom yes rot yes loop all
#Setup output
thermo_style custom step temp press etotal vol
thermo 100
#Colored thermal bath
fix scapegoat_qtb all nve #NVE does the time integration
fix methane_qtb all qtb temp ${temperature} damp ${damp_qtb} seed 35082 f_max 0.3 N_f 50 #Change f_max if your Debye frequency is higher
timestep ${delta_t}
run 3000 #750 fs
unfix methane_qtb
unfix scapegoat_qtb