Release of "learned_free_energy_estimation".

PiperOrigin-RevId: 336058713
2026-05-09 21:07:49 +08:00 · 2020-10-08 12:28:23 +01:00
parent c82e368e0e
commit 1d86410c90
7 changed files with 900936 additions and 0 deletions
@@ -0,0 +1,129 @@
+# 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}