Hi all,
Following is my code I use to run NEMD after NVE. Since the NEMD process is long and there is a 3 day time limit on the GPU, I have to restart the simulation. However, I am encountering a problem where if I restart where I left off using the exact same script but different timestep, the fluid in the domain would become frozen. Is there a solution to this or does this have to do with the fixes?
---------- Initialize Simulation ----------------------------------------------
read_restart ${restart_file}
#-------------------Define Simulation Box -----------------------------------------------
region box block 0 2508.8 0 ${domain_ymax} 0 39.201 units box
#-------------------Define Parameters -------------------------------------------
variable Sig_Ar equal 3.405 # in Angstrom
variable cutoff equal {NSigma}*{Sig_Ar}
#-------------------Define (Initial) Argon Region -----------------------------------
mass 1 39.948002 #Ar
#-------------------Wall Regions -----------------------------------------------
lattice fcc 3.92 #Pt
region lowrwall block 156.30 470.90 0.00 {lowerwall_ymax} -0.48 38.72 units box
region uperwall block 156.30 470.90 {upperwall_ymin} ${upperwall_ymax} -0.48 38.72 units box
region lowrwallwr block 156.30 470.90 0.00 1.96 0.00 39.201 units box
region uperwallupr block 156.30 470.90 {upperfr_ymin} {upperfr_ymax} 0.00 39.201 units box # frozen region
region hotlwr block 156.30 233.74 {lowerwall_thermo_ymin} {lowerwall_thermo_ymax} -0.48 38.72 units box
region hotupr block 156.30 233.74 {upperwall_thermo_ymin} {upperwall_thermo_ymax} -0.48 38.72 units box
region coldlwr block 393.46 470.90 {lowerwall_thermo_ymin} {lowerwall_thermo_ymax} -0.48 38.72 units box
region coldupr block 393.46 470.90 {upperwall_thermo_ymin} {upperwall_thermo_ymax} -0.48 38.72 units box
region hotlwrnve block 156.30 233.74 {lowerwall_thermo_ymax} {lowerwall_ymax} -0.48 38.72 units box
region hotuprnve block 156.30 233.74 {upperwall_ymin} {upperwall_thermo_ymin} -0.48 38.72 units box
region coldlwrnve block 393.46 470.90 {lowerwall_thermo_ymax} {lowerwall_ymax} -0.48 38.72 units box
region colduprnve block 393.46 470.90 {upperwall_ymin} {upperwall_thermo_ymin} -0.48 38.72 units box
-----------------Define XYZ Region for Particle Tracking ---------------------------
region Track block 480.00 600.00 0 92.121 0 39.201 units box
mass 2 195.09 #Pt
#------------------Divide Groups------------------------------------------
group argon type 1
group pt type 2
group glowrwall region lowrwall
group guperwall region uperwall
group ptfrzlwr region lowrwallwr
group ptfrzupr region uperwallupr
group ghotlwr region hotlwr
group ghotupr region hotupr
group ghotlwrnve region hotlwrnve
group ghotuprnve region hotuprnve
group gcoldlwr region coldlwr
group gcoldupr region coldupr
group gcoldlwrnve region coldlwrnve
group gcolduprnve region colduprnve
group hotwall union ghotlwr ghotupr
#group hotwallnve union ghotlwrnve ghotuprnve
group coldwall union gcoldlwr gcoldupr
#group coldwallnve union gcoldlwrnve gcolduprnve
group applynve union argon ghotlwrnve ghotuprnve gcoldlwrnve gcolduprnve
---------- Define Interatomic Potential------------------------------------
pair_style hybrid lj/cut {cutoff} eam
pair_coeff 1 1 lj/cut 0.0104232 3.405 #Cao,Chen,Guo
pair_coeff 2 2 eam Pt_u3.eam
variable applied_eps equal 0.0055798*{eps}
pair_coeff 1 2 lj/cut ${applied_eps} 3.085 #Cao,Chen,Guo
neighbor 4 bin #for different types of atoms? 4 means 4sigma
neigh_modify delay 10 check no
---------- Computes for PE, KE, T & Stress ------------------------------------
compute ar_temp argon temp
compute hotwall_temp hotwall temp
compute coldwall_temp coldwall temp
compute pe argon pe/atom
compute ke argon ke/atom
---------- Heat Flux Calculations ----------------------------------------------------------------------------
variable energy atom c_ke+c_pe
--------------------------------------------------------------------------------------------------------------
compute heat_compute all chunk/atom bin/2d y 0.98 0.392 x lower 0.392 bound y 1.08 51.84 bound x 0.100 313.5 units box #0_1_sigma (to prevent the inclusion of extra bins, span is shortened/streched by 0.1A)
--------------------------------------------------------------------------------------
#fix 1 all nve
fix 93 applynve nve
fix 94 hotwallnve nve
fix 95 coldwallnve nve
#variable applied_heat equal {nW}*0.0062415# fix 101 hotwall nvt temp {temperature} {THot} 0.01 fix 102 coldwall nvt temp {temperature} ${TCold} 0.01
compute 11 argon chunk/atom bin/2d y lower 0.392 x lower 0.392 bound y 1.08 91.04 bound x 0.100 627.1 units box
compute 11 argon chunk/atom bin/2d y lower 0.392 x lower 0.392 bound y 0 ${domain_ymax} bound x 0 627.2 units box #0_1_sigma (to prevent the inclusion of extra bins, span is shortened/streched by 0.1A)
#fix energyRampT argon ave/chunk 4 100000 400000 11 c_pe c_ke file energy_rampT.profile
thermo_style custom step cpu c_ar_temp c_hotwall_temp c_coldwall_temp
thermo 10000
timestep 0.005 # dump 2 all xyz 50000 NEMD.xyz
dump 2 all custom ${freq} coords_transient.dump id type x y z
restart 1000000 NEMD.restart
run 10000000
#------------------- unfix from transient
unfix 101
unfix 102
-------------------steady state simulations
fix 103 hotwall nvt temp {THot} {THot} 0.01
fix 104 coldwall nvt temp {TCold} {TCold} 0.01
---------- Computes for PE, KE, T & Stress ------------------------------------
compute temp all temp
compute stress all stress/atom temp
compute stressv all stress/atom NULL virial
compute 11 argon chunk/atom bin/2d y lower 0.392 x lower 0.392 bound y 0.48 106.32 bound x 0.100 627.1 units box
---------- fixes ----------------------------------------------------------------------------
fix 110 argon ave/chunk 1 400000 400000 11 density/mass file dens_nemd_01_sig.profile
fix 111 argon ave/chunk 1 400000 400000 11 c_pe c_ke file energy_nemd_01_sig.profile
fix 113 argon ave/chunk 1 400000 400000 11 c_stress[] file stress_nemd_01_sig.profile
fix 114 argon ave/chunk 1 400000 400000 11 c_stressv[] file stressv_nemd_01_sig.profile
fix 115 argon ave/chunk 1 400000 400000 11 vx vy vz file vel_nemd_01_sig.profile
fix 1101 argon ave/chunk 1 400000 400000 11 temp file temp_nemd_01_sig.profile
dump 9 all custom ${freq} coords.dump id type x y z
restart 1000000 NEMD.restart
run 24000000