#NEMD simulation, Constant Heat Flux #NOTE: For a description of any command, Google "LAMMPS command". #------------SET SIM SEED-------------------------------------------------- variable seed equal 64576 #------------READ STRUCTURE------------------------------------------------ units metal #Use Metal style units atom_style full #Use 'molecular' style for solids,molecules boundary p p p #Periodic BC on x, y, z read_data sam_coord_in.data #Use lattice structure from file group lfix molecule 1 #Left fixed atoms group hotres molecule 2 #Host reservoir atoms group coldres molecule 29 #Cold reservoir atoms group rfix molecule 30 #Right fixed atoms group sample subtract all lfix hotres coldres rfix #Sample atoms group moving union hotres sample coldres #All non-fixed atoms group fixed molecule 1 30 #All fixed atoms group temp_nemd type 1 4 group shake_water type 5 6 group ox type 5 group hy type 6 set group ox charge -1.0484 set group hy charge 0.5242 #------------POTENTIALS------------------------------------------- pair_style hybrid/overlay eam lj/cut/coul/long 10.0 8.5 lj/cut/tip4p/long 5 6 3 3 0.125 10.0 8.5 morse 9 special_bonds lj 0.0 0.0 0.0 pair_coeff 1 1 eam /home/shubhadm/lammps-17Jun13/potentials/Au_u3.eam #EAM potential interaction for Au-Au pair_coeff 4 4 eam /home/shubhadm/lammps-17Jun13/potentials/Au_u3.eam #EAM potential interaction for Au-Au pair_coeff * * morse ${M_Zero} 1.583000 3.024200 #Morse pair potential interaction for all pair_coeff 1 5 lj/cut/coul/long 0.02558 3.6000 #LJ pair potential interaction for Au-O pair_coeff 4 5 lj/cut/coul/long 0.02558 3.6000 #LJ pair potential interaction for Au-O pair_coeff 5 5 lj/cut/tip4p/long 0.0070575 3.16435 #LJ pair potential interaction for O-O pair_coeff 5 6 lj/cut/tip4p/long 0.0 0.0 #LJ pair potential interaction for O-H pair_coeff 6 6 lj/cut/tip4p/long 0.0 0.0 #LJ pair potential interaction for H-H pair_coeff 1 6 lj/cut/coul/long 0.0 0.0 #LJ pair potential interaction for Au-H pair_coeff 4 6 lj/cut/coul/long 0.0 0.0 #LJ pair potential interaction for Au-H pair_coeff 1 2 morse 0.380000 1.470000 2.650000 #Morse pair potential interaction for S-Au pair_coeff 2 4 morse 0.380000 1.470000 2.650000 #Morse pair potential interaction for S-Au pair_coeff 2 2 lj/cut/coul/long 0.017240 4.250000 #LJ pair potential interaction for S-S pair_coeff 3 3 lj/cut/coul/long 0.005130 3.914000 #LJ pair potential interaction for C-C pair_coeff 1 3 lj/cut/coul/long 0.002900 3.424000 #LJ pair potential interaction for Au-C pair_coeff 3 4 lj/cut/coul/long 0.002900 3.424000 #LJ pair potential interaction for Au-C pair_coeff 2 3 lj/cut/coul/long 0.007500 3.732000 #LJ pair potential interaction for C-S pair_coeff 2 5 lj/cut/coul/long 0.008800 3.357200 #LJ pair potential interaction for S-O pair_coeff 2 6 lj/cut/coul/long 0.0 0.0 #LJ pair potential interaction for S-H pair_coeff 3 5 lj/cut/coul/long 0.00600 3.539200 #LJ pair potential interaction for C-O pair_coeff 3 6 lj/cut/coul/long 0.0 0.0 #LJ pair potential interaction for C-H pair_modify shift yes #------------DEFINE HARMONIC POTENTIAL FOR C-C BOND------------------------------------------- bond_style harmonic bond_coeff 1 11.27 1.54 #------------DEFINE HARMONIC POTENTIAL FOR C-S BOND------------------------------------------- bond_coeff 2 9.54 1.815 #------------DEFINE HARMONIC POTENTIAL FOR O-H BOND------------------------------------------- bond_coeff 3 19.51 0.9572 #------------DEFINE HARMONIC ANGLE POTENTIAL FOR C-C-C BOND------------------------------------------- angle_style harmonic angle_coeff 1 2.6940 109.5 #------------DEFINE HARMONIC ANGLE POTENTIAL FOR S-C-C BOND------------------------------------------- angle_coeff 2 2.6940 114.4 #------------DEFINE HARMONIC ANGLE POTENTIAL FOR H-O-H BOND------------------------------------------- angle_coeff 3 2.39 104.5 #------------DEFINE TORSION ANGLE POTENTIAL FOR S-C-C-C and C-C-C-C BONDS----------------------------- dihedral_style multi/harmonic dihedral_coeff 1 0.09617 -0.125988 -0.13598 0.0317 0.27196 0.32642 kspace_style pppm 1.0e-5 #final npt relaxation neighbor 2.0 bin neigh_modify delay 0 every 1 check yes #------------DEFINE VARIABLES---------------------------------------------- #------------thermo Parameters variable T_run equal 300 #Run temperature (K) variable heat equal 0.1 #The heat to be added and removed per time step (eV/ps) variable dt equal 5e-4 #The simulation timestep (ps) #--------DEFINE CUSTOM COMPUTES---------------------------------------- compute movingTemp moving temp #Define a compute to compute the temperature of moving atoms compute stressPerAtom sample stress/atom #Define a compute to compute the stress per atom of sample atoms compute p sample reduce sum c_stressPerAtom[1] c_stressPerAtom[2] c_stressPerAtom[3] #Compute the components of pressure of sample atoms log log_300K_NVT.lammps #--------NVT RESCALE -------------------------------------------------- velocity moving create ${T_run} ${seed} rot yes dist gaussian #Set initial velocities fix 2 moving nvt temp ${T_run} ${T_run} 10 #Fix NVT ensemble fix 1 shake_water shake 1e-6 200 100000 b 3 a 3 variable samplePressX equal -(c_p[1])/(vol*7/10) variable samplePressY equal -(c_p[2])/(vol*7/10) variable samplePressZ equal -(c_p[3])/(vol*7/10) thermo_style custom step c_movingTemp v_samplePressX v_samplePressY v_samplePressZ press etotal pe ke vol thermo 1000 timestep ${dt} run 500000 unfix 2 #Unfix NVT #--------CONSTANT HEAT FLUX NEMD--------------------------------------- fix hot_res hotres heat 1 ${heat} #At every time step, add energy 'heat' to hot reservoir fix cold_res coldres heat 1 -${heat} #At every time step, remove energy 'heat' from cold reservoir fix 2 all nve #Fix NVE ensemble fix lwall lfix setforce 0.0 0.0 0.0 #Set force on left fixed atoms to 0 fix rwall rfix setforce 0.0 0.0 0.0 #Set force on right fixed atoms to 0 variable samplePressNet2 equal -(c_p[1]+c_p[2]+c_p[3])/(3*vol*7/10) #Compute net pressure using custom computation thermo_style custom step c_movingTemp v_samplePressNet2 press etotal pe ke vol thermo 1000 timestep ${dt} run 2500000 log log_au_sam_water_300K_1.lammps #--------CONSTANT HEAT FLUX NEMD--------------------------------------- compute 3 temp_nemd ke/atom compute 4 temp_nemd atom/molecule c_3 fix avetemp1 all ave/time 100 1 100 & c_4[1] c_4[2] c_4[3] c_4[4] c_4[5] c_4[6] c_4[7] c_4[8] c_4[9] c_4[10] c_4[11] & c_4[12] c_4[13] c_4[14] c_4[15] c_4[16] c_4[17] c_4[18] & file log_temp_stand_300K_10c2s_eam2_large_new_cut9_water_tip4plong_4.lammps variable samplePressNet2 equal -(c_p[1]+c_p[2]+c_p[3])/(3*vol*7/10) #Compute net pressure using custom computation thermo_style custom step c_movingTemp v_samplePressNet2 press etotal pe ke vol #Log thermo 1000 dump 5 all xyz 1000000 sam_stand_300K_10c2s_min_large_fullstruct_cut9_water_4.xyz timestep ${dt} run 20000000 unfix 2 unfix 1 #Unfix NVE unfix cold_res #Unfix cold energy removal unfix hot_res #Unfix hot energy addition unfix lwall #Unfix left atoms unfix rwall #Unfix right atoms uncompute movingTemp uncompute stressPerAtom uncompute p uncompute 3 uncompute 4