# sbatch --account=XXXXX --mail-user=XXXXX --mail-type=ALL --ntasks=1 --mem=4G --time=1:00:00 --array="0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0" --wrap "srun lmp_beluga -in calf15.in -var pressure_bar $SLURM_ARRAY_TASK_ID" # The state point # temperature 273.15 variable temp index 273.15 # K variable press index 1.0 # bar units real atom_style full boundary p p p # specify cutoffs # cutoff 6.39 angstrom # delr 1.0 angstrom ??? MAXIMUM MC MOVE ??? #pair_style lj/cut/coul/long 6.39 pair_style lj/cut 6.39 # dealing with coulombic forces # ewald precision 1d-6 #kspace_style ewald 1e-6 lattice custom 1.0 a1 13.8034 0 0 a2 0 12.9639 0 a3 0 0 16.8318 basis 0 0 0 read_data calf15.data molecule co2mol CO2.txt group co2 type 3 9 neighbor 2.0 bin neigh_modify every 1 delay 10 check yes velocity all create ${temp} 54654 timestep 1.0 pair_coeff Nc Nc 0.06899 3.260 pair_coeff Nc Nh 0.06899 3.260 pair_coeff Nc Cc 0.08511 3.345 pair_coeff Nc C 0.08511 3.345 pair_coeff Nc Co 0.064056 3.0025 pair_coeff Nc Oc 0.108369 3.1385 pair_coeff Nc H5 0.05509 2.915 pair_coeff Nc Hn 0.05509 2.915 pair_coeff Nc O 0.06433 3.190 pair_coeff Nc Zn 0.09249 2.860 pair_coeff Nh Nh 0.06899 3.260 pair_coeff Nh Cc 0.08511 3.345 pair_coeff Nh C 0.08511 3.345 pair_coeff Nh Co 0.064056 3.0025 pair_coeff Nh Oc 0.108369 3.1385 pair_coeff Nh H5 0.05509 2.915 pair_coeff Nh Hn 0.05509 2.915 pair_coeff Nh O 0.06433 3.190 pair_coeff Nh Zn 0.09249 2.860 pair_coeff Cc Cc 0.10499 3.430 pair_coeff Cc C 0.10499 3.430 pair_coeff Cc Co 0.079023 3.0875 pair_coeff Cc Oc 0.133689 3.2235 pair_coeff Cc H5 0.06796 3.000 pair_coeff Cc Hn 0.06796 3.000 pair_coeff Cc O 0.07936 3.275 pair_coeff Cc Zn 0.11410 2.945 pair_coeff C C 0.10499 3.430 pair_coeff C Co 0.079023 3.0875 pair_coeff C Oc 0.133689 3.2235 pair_coeff C H5 0.06796 3.000 pair_coeff C Hn 0.06796 3.000 pair_coeff C O 0.07936 3.275 pair_coeff C Zn 0.11410 2.945 pair_coeff Co Co 0.059477 2.7450 pair_coeff Co H5 0.051152 2.6575 pair_coeff Co Hn 0.051152 2.6575 pair_coeff Co O 0.059732 2.9325 pair_coeff Co Oc 0.100621 2.8810 pair_coeff Co Zn 0.085875 2.6025 pair_coeff H5 H5 0.04399 2.570 pair_coeff H5 Hn 0.04399 2.570 pair_coeff H5 O 0.05137 2.845 pair_coeff H5 Oc 0.086537 2.7935 pair_coeff H5 Zn 0.07386 2.515 pair_coeff Hn Hn 0.04399 2.570 pair_coeff Hn O 0.05137 2.845 pair_coeff Hn Oc 0.086537 2.7935 pair_coeff Hn Zn 0.07386 2.515 pair_coeff O O 0.05999 3.120 pair_coeff O Oc 0.101052 3.0685 pair_coeff O Zn 0.08624 2.790 pair_coeff Oc Oc 0.170228 3.0170 pair_coeff Oc Zn 0.145280 2.7385 pair_coeff Zn Zn 0.12399 2.460 ############################### FIX GCMC ############################### # fix ID group-ID gcmc N X M type seed T mu displace keyword values ... # ID, group-ID are documented in fix command # gcmc = style name of this fix command # N = invoke this fix every N steps # X = average number of GCMC exchanges to attempt every N steps # M = average number of MC moves to attempt every N steps # type = atom type for inserted atoms (must be 0 if mol keyword used) # seed = random # seed (positive integer) # T = temperature of the ideal gas reservoir (temperature units) # mu = chemical potential of the ideal gas reservoir (energy units) # displace = maximum Monte Carlo translation distance (length units) # zero or more keyword/value pairs may be appended to args # mu (CO2) = -18 KJ/mol = -4.30 kcal/mol # Rahbari, A., Poursaeidesfahani, A., Torres-Knoop, A., Dubbeldam, D., & Vlugt, T. J. H. # (2018). Chemical potentials of water, methanol, carbon dioxide and hydrogen sulphide at low # temperatures using continuous fractional component Gibbs ensemble Monte Carlo. Molecular # Simulation, 44(5), 405-414. https://doi.org/10.1080/08927022.2017.1391385] # As an alternative to specifying mu directly, the ideal gas reservoir can be # defined by its pressure P using the pressure keyword, in which case the # user-specified chemical potential is ignored. # *** LAMMPS real units pressure atm, 1 bar = 1.013 atm *** # rigid constraints with thermostat fix myrigid co2 rigid/small molecule mol co2mol # dynamically update fix rigid/nvt/small temperature ndof fix_modify myrigid dynamic/dof yes variable tfac equal 5.0/3.0 # (3 trans + 2 rot)/(3 trans) fix 1 co2 gcmc 1 1 0 0 12345 ${temp} 0.0 1.0 & pressure $(v_press/1.013) & mol co2mol & tfac_insert ${tfac} & group co2 rigid myrigid # Use of this fix typically will cause the number of atoms to fluctuate, therefore, you # will want to use the compute_modify dynamic/dof command to ensure that the current # number of atoms is used as a normalizing factor each time temperature is computed. ############################### FIX SETFORCE ############################### # freeze all non-CO2 atoms group mof subtract all co2 fix 2 mof setforce 0.0 0.0 0.0 ############################### COMPUTE TEMP ############################### # If used with fix nvt, the temperature of the imaginary reservoir, T, # should be set to be equivalent to the target temperature used in fix nvt. # Otherwise, the imaginary reservoir will not be in thermal equilibrium with # the simulation cell. Also, it is important that the temperature used by # fix nvt is dynamically updated, which can be achieved as follows: compute mdtemp all temp compute_modify mdtemp dynamic/dof yes fix mdnvt all nvt temp ${temp} ${temp} 10.0 fix_modify mdnvt temp mdtemp ######################################################################## # determine how often to write to the 'numguests.out' file # numguests 100 variable numguests equal count(co2) thermo_style custom step temp press pe ke density atoms v_numguests thermo 10 # unable to turn off title line so use empty string # https://matsci.org/t/fix-print-without-title/25204 fix 3 all print 100 "$(step) ${numguests}" file calf15-numguests.out screen no title "" # store GCMC frames for visualization every [n] steps # history 10000 dump 1 all custom 1 calf15.lammpstrj id mol element type x y z ix iy iz ############################## EQUILIBRATION ############################## # Number of MC steps to perform to equilibrate the system. # After this number of MC steps, averaging occurs. #run 500000 run 20000 ############################### PRODUCTION ############################### # Number of steps to perform after equilibration. # The simulation will stop after this number of steps has been reached. #run 500000 run 20000