Hello All,
I am simulating the evaporation of water in a capacitor by utilizing a 2D simulation. Before I institute any electric field constraints I attempt to equilibrate the system to a temperature of 300K. I use a velocity create command to make the temperature of the capacitor walls (which are gold) and water 300K (I constrain the water using fix RATTLE with bond and angle constraints creating 1000 frozen angles), utilize a run 0 command and then utilize a velocity scale command to explicitly make the temperature of the capacitor walls and water molecules 300K. I utilize a time step of 0.1 femtoseconds. I apply NVT constraints to the capacitor walls and water for 500 picoseconds (5,000,000 time steps) and then I apply NVE constraints to the capacitor walls and water for another 500ps. I compute the temperature of the capacitor walls and water (separately) every 100 time steps. The capacitor and water temperature oscillate about a mean of 300K for the first 500ps but upon initiation of the NVE constraints the water temperature rises to around 350K while the capacitor walls are still oscillaing about that mean temperature of 300K. I would like to know why the temperature of my water molecules just rises when I switch from an NVT fix to an NVE fix? Any Insight will be greatly appreciated and I have attached part of my input script for reference. Thank you.
Sincerely,
Malcolm Porterfield
###Overall Simulation Parameters###
clear #Deletes all atoms, restores all settings to their default values, and frees all memory allocated by LAMMPS
echo both #Output code to both log file and terminal screen
units metal #Units are of metal style (distance=Angstrom, time=picosecond, energy=eV, temperature=K, pressure=bars)
atom_style full #Atoms have charge as well a molecular attributes (i.e bonds, angles, dihedrals and impropers)
boundary p f p #non-periodic and fixed boundary condition in the y-direction, periodic boundary conditions in the x and z directions
dimension 2 #Create a 2-dimensional simulation
lattice sq2 4.08 #Create a lattice constant of 4.08 Angstroms sq2 stlye (because the simulation is two dimensional) style to be applied to Gold
#Gold has an FCC style lattice with 4 basis atoms (1 at the corner and 3 at the cube face centers). The sq2 style has 2 basis
#atoms in total, 1 at the corner of the square and 1 atom in the center
###Setting Atomic interactions/Attributes###
##Utilizing Simple Point Charge water molecules as given in LAMMPS manual##
mass 1 1.008 #Set mass of type 1 atoms (Hydrogen)
mass 2 15.9994 #Set mass of type 2 atoms (Oxygen)
mass 3 196.7 #Set the mass of type 3 atoms (Gold atoms in the left side of Capacitor)
mass 4 196.7 #Set the mass of type 4 atoms (Gold atoms in the right side of Capacitor)
mass 5 196.7 #Set the mass of type 5 atoms (Gold atoms on the outer edge of each side of the Capacitor)
set type 3*5 charge 0.0 #Set the charge of type 3 through type 5 atoms to zero
pair_style hybrid/overlay eam/alloy lj/cut 9 coul/cut 9 morse 11 #Allow for multiple pair styles to be used in general and per interaction
#pair
#Pair style eam/alloy will be used for the Gold to Gold atom interactions (embedded atom method)
#Pair style lj/cut uses a standard lennard-jones interaction potential. The cutoff is 9 angstroms.
#Pair style coul/cut uses a standard Coulombic pairwise interaction given by Coulomb's Law. The cutoff is 9 angstroms.
#Pair style morse uses a standard morse potential interaction with a cutoff value of 11 angstroms
#Gold
pair_coeff * * eam/alloy Au.eam.alloy NULL NULL Au Au Au #Maps the atom types 3, 4 and 5 to the Au element in the setfl file
#Water
pair_coeff 1 2 coul/cut #hydrogen-oxygen
pair_coeff 1 1 coul/cut #hydrogen-hydrogen
pair_coeff 2 2 coul/cut #oxygen-oxygen part 1
pair_coeff 2 2 lj/cut 0.006734 3.166 #oxygen-oxygen part 2
#Gold-water
pair_coeff 1 3*5 morse 0.00082914504 1.41 4.14 #Interaction between hydrogen and gold (D0, alpha, r0 in that order for the
#constants)
pair_coeff 2 3*5 morse 0.01927754 0.905 4.2 #Interaction between oxygen and gold
#Bond Parameters
bond_style harmonic #Harmonic bond style
bond_coeff 1 100 1.0 #Apply to type 1 bonds with a K factor of 100 eV/m^2 (I believe this makes the bonds more stiff) and equilibrium
#length of 1 Angstrom
#Angle Parameters
angle_style harmonic #Harmonic angle style
angle_coeff 1 100 109.47 #Apply to type 1 angles with a K factor 100 rad/m^2 (I believe this makes the bonds more stiff) and
#equilibrium angle of 109.47 degrees
###Simulation Initializations###
velocity LeftCapg create 300 123456 #Effect the kinetic temperature of the atoms in the LeftCapg group
#set the temperature to 300K initially (123456 is a random seed number)
velocity RightCapg create 300 123457 #Effect the kinetic temperature of the atoms in the RightCapg group
#set the temperature to 300K initially (123456 is a random seed number)
velocity Water create 300 124568 #Effect the kinetic temperature of the atoms in the Water group
#set the temperature to 300K initially (124568 is a random seed number)
fix constrain Water rattle 1.0e-4 100 1000 b 1 a 1 #Create a fix called constrain that effects the atoms in the water group that, employs
#the SHAKE algorithm
#has a solution tolerance of 1e-4, has a value of 100 as the maximum number of iterations
#to find a solution,
#Prints the SHAKE statistics every 1000 time steps, utilizes one bond type and one angle
#type
run 0 #ensuring all degrees of freedom are properly acconted for
velocity LeftCapg scale 300 #Explicitly scale the temperature of the LeftCapg atoms to 300K
velocity RightCapg scale 300 #Explicitly scale the temperature of the RighCapg atoms to 300K
velocity Water scale 300 #Explicitly scale the temperature of the Water molecules to 300K
timestep 0.0001 #Have a simulation time step of 0.0001 picoseconds (or 0.1 femtosecond)
neighbor 2.0 bin #Have a skin distance of 2 Angstroms and sort the atoms by binning
fix Pool all wall/reflect ylo EDGE #Create a fix called Pool that effects all atoms and will reflect atoms that attempt to move through
#the bottom of the simulation domian (bottom y-axis)
fix Pool2 all wall/reflect yhi EDGE #Create a fix called Pool2 that effects all atoms and will reflect atoms that attempt to move
#through
#the bottom of the simulation domian (top y-axis)
###System Temperature Equilibration###
unfix constrain
fix 1 Water nvt temp 300 300 0.1 #create a fix called 1 that effects the atoms in the Water group with a canonical ensemble
#thermostat
#have the temperature start at 300K and end at 300K with a damping parameter of 0.1 picosecond
fix 2 LeftCapg nvt temp 300 300 0.1 #create a fix called 2 that effects the atoms in the LeftCapg group with a canonical
#ensemble thermostat
#have the temperature start at 300K and end at 300K with a damping parameter of 0.1
#picosecond
fix 3 RightCapg nvt temp 300 300 0.1 #create a fix called 3 that effects the atoms in the RightCapg group with a canonical
#ensemble thermostat
#have the temperature start at 300K and end at 300K with a damping parameter of 0.1
#picosecond
fix 2dsim all enforce2d #create a fix called 2dsimbef that effects all atoms and zeros out the z-component of their velocity and force
fix constrain Water rattle 1.0e-4 100 1000 b 1 a 1 #Create a fix called constrain that effects the atoms in the water group that,
#employs the SHAKE algorithm
#has a solution tolerance of 1e-4, has a value of 100 as the maximum number of
#iterations to
#find a solution,
#Prints the SHAKE statistics every 1000 time steps, utilizes one bond type and one
#angle
#type
run 5000000