forked from lijiext/lammps
git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@13682 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
sjplimp
parent
92b5c350d8
commit
6864dfa2d1
|
|
@ -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.
|
||||
|
|
@ -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
|
||||
|
|
@ -0,0 +1,55 @@
|
|||
## 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}
|
||||
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}
|
||||
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}
|
||||
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 10000 #10 ps
|
||||
|
|
@ -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} #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 100
|
||||
run 2000 # 2 ps
|
||||
unfix quartz_qtb
|
||||
unfix scapegoat_qtb
|
||||
File diff suppressed because it is too large
Load Diff
|
|
@ -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} #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 100
|
||||
run 20000 # 20 ps
|
||||
unfix quartz_qtb
|
||||
unfix scapegoat_qtb
|
||||
|
|
@ -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
|
||||
|
|
@ -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
|
||||
|
|
@ -0,0 +1,63 @@
|
|||
## 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 10.0 1 1 1.0e-6
|
||||
pair_coeff * * ffield.reax 1 2
|
||||
|
||||
#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
|
||||
#Initilization
|
||||
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
|
||||
|
|
@ -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
|
||||
|
|
@ -0,0 +1,69 @@
|
|||
## 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 10.0 1 1 1.0e-6
|
||||
pair_coeff * * ffield.reax 1 2
|
||||
|
||||
#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
|
||||
#Initilization
|
||||
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
|
||||
Loading…
Reference in New Issue