mirror of
https://github.com/google-deepmind/deepmind-research.git
synced 2026-05-09 21:07:49 +08:00
Peter Wirnsberger
130 lines
5.3 KiB
Plaintext
130 lines
5.3 KiB
Plaintext
# Simulate a WCA solute particle immersed into a Lennard-Jones solvent.
|
|
#
|
|
# We first generate an initial configuration of particles located at random
|
|
# positions. We then minimise the system and equilibrate it to the target
|
|
# temperature. During the production run, we sample configurations at regular
|
|
# intervals.
|
|
|
|
################################################################################
|
|
# Specify input parameters.
|
|
################################################################################
|
|
|
|
# We use the same box dimensions and temperature as in Jarzynski (2002),
|
|
# Phys. Rev. E 65, 046122, but convert these quantities to reduced units.
|
|
# T=300 K, eps=0.1854 kcal/mol, k_boltzmann=0.001987204 kcal/mol/K.
|
|
variable solute_radius equal 9.2/3.542
|
|
variable box_length equal 22.28/3.542
|
|
variable temperature_reference equal 300*0.001987204/0.1854
|
|
variable temperature_damp equal 0.5
|
|
variable seed equal 1234
|
|
variable num_solvent_particles equal 125
|
|
variable dump_frequency equal 5000
|
|
variable timestep equal 0.002
|
|
variable runtime_equi equal 50000
|
|
variable runtime_prod equal 500000
|
|
variable steps_equilibration equal floor(${runtime_equi}/${timestep})
|
|
variable steps_production equal floor(${runtime_prod}/${timestep})
|
|
|
|
################################################################################
|
|
# Create solvent particles at random positions and place solute at the origin.
|
|
################################################################################
|
|
|
|
units lj
|
|
atom_style atomic
|
|
dimension 3
|
|
boundary p p p
|
|
variable box_length_half equal ${box_length}/2.0
|
|
region domain block -${box_length_half} ${box_length_half} &
|
|
-${box_length_half} ${box_length_half} &
|
|
-${box_length_half} ${box_length_half} &
|
|
units box
|
|
create_box 2 domain
|
|
create_atoms 2 random 1 ${seed} domain
|
|
create_atoms 1 random ${num_solvent_particles} ${seed} domain
|
|
set type 2 x 0.0 y 0.0 z 0.0
|
|
group solute type 2
|
|
group solvent type 1
|
|
|
|
# Initialise particle masses.
|
|
mass 1 1
|
|
mass 2 10000000
|
|
|
|
################################################################################
|
|
# Define pair-style.
|
|
################################################################################
|
|
|
|
variable sigma_solute equal ${solute_radius}
|
|
variable radius_cutoff_solute equal 1.12246*${sigma_solute}
|
|
variable radius_cutoff equal ${box_length}/2.0
|
|
pair_style lj/cut ${radius_cutoff}
|
|
|
|
# Lennard-Jones interaction between solvent particles.
|
|
pair_coeff 1 1 1.0 1.0
|
|
|
|
# WCA interaction between solute and solvent particles.
|
|
pair_coeff 1 2 1.0 ${sigma_solute} ${radius_cutoff_solute}
|
|
|
|
# No interaction between solute particles.
|
|
pair_coeff 2 2 0.0 0.0
|
|
|
|
# Shift energies to be zero at the cutoff and update neighbor list settings.
|
|
pair_modify shift yes
|
|
neighbor 0.3 bin
|
|
neigh_modify delay 5
|
|
|
|
################################################################################
|
|
# Perform minimisation.
|
|
################################################################################
|
|
|
|
fix freeze solute setforce 0.0 0.0 0.0
|
|
minimize 1.0e-4 1.0e-6 100 1000
|
|
|
|
################################################################################
|
|
# Perform equilibration run.
|
|
################################################################################
|
|
|
|
reset_timestep 0
|
|
|
|
compute temperature_solvent solvent temp
|
|
compute kinetic_energy all ke
|
|
compute potential_energy all pe
|
|
variable total_energy equal c_kinetic_energy+c_potential_energy
|
|
|
|
# Compute centre-of-mass velocity and solute position for monitoring.
|
|
variable vcmx equal "vcm(all,x)"
|
|
variable vcmy equal "vcm(all,y)"
|
|
variable vcmz equal "vcm(all,z)"
|
|
variable vcm2 equal v_vcmx*v_vcmx+v_vcmy*v_vcmy+v_vcmz*v_vcmz
|
|
compute position_solute solute com
|
|
|
|
# Specify terminal output.
|
|
thermo_style custom step temp c_temperature_solvent c_potential_energy &
|
|
c_kinetic_energy v_total_energy press v_vcm2 &
|
|
c_position_solute[1] c_position_solute[2] c_position_solute[3]
|
|
thermo_modify norm no
|
|
thermo 1000
|
|
|
|
# Specify integration timestep.
|
|
timestep ${timestep}
|
|
|
|
# Apply Langevin thermostat to solvent particles.
|
|
fix flangevin solvent langevin ${temperature_reference} &
|
|
${temperature_reference} ${temperature_damp} ${seed} zero yes
|
|
|
|
fix fnve all nve
|
|
run ${steps_equilibration}
|
|
|
|
################################################################################
|
|
# Production run.
|
|
################################################################################
|
|
|
|
reset_timestep 0
|
|
|
|
dump fDumpTrajectory all custom ${dump_frequency} &
|
|
trajectory.dat id type x y z
|
|
dump_modify fDumpTrajectory sort id format float %.8g
|
|
|
|
fix fDumpEnergy all ave/time ${dump_frequency} 1 ${dump_frequency} &
|
|
c_potential_energy file energy.dat format %.8g
|
|
run ${steps_production}
|