|
| 1 | +from ase.units import kJ, mol |
| 2 | +import numpy as np |
| 3 | + |
| 4 | +import psiflow |
| 5 | +from psiflow.data import Dataset |
| 6 | +from psiflow.geometry import Geometry |
| 7 | +from psiflow.hamiltonians import PlumedHamiltonian, MACEHamiltonian |
| 8 | +from psiflow.sampling import Walker, sample, quench, Metadynamics, replica_exchange |
| 9 | + |
| 10 | + |
| 11 | +PLUMED_INPUT = """UNITS LENGTH=A ENERGY=kj/mol |
| 12 | +d_C: DISTANCE ATOMS=3,5 |
| 13 | +d_O: DISTANCE ATOMS=1,5 |
| 14 | +CV: COMBINE ARG=d_C,d_O COEFFICIENTS=1,-1 PERIODIC=NO |
| 15 | +
|
| 16 | +""" |
| 17 | + |
| 18 | + |
| 19 | +def get_bias(kappa: float, center: float): |
| 20 | + plumed_str = PLUMED_INPUT |
| 21 | + plumed_str += '\n' |
| 22 | + plumed_str += 'RESTRAINT ARG=CV KAPPA={} AT={}\n'.format(kappa, center) |
| 23 | + return PlumedHamiltonian(plumed_str) |
| 24 | + |
| 25 | + |
| 26 | +def main(): |
| 27 | + aldehyd = Geometry.load('data/acetaldehyde.xyz') |
| 28 | + alcohol = Geometry.load('data/vinyl_alcohol.xyz') |
| 29 | + |
| 30 | + mace = MACEHamiltonian.mace_cc() |
| 31 | + energy = mace.compute([aldehyd, alcohol], 'energy').result() |
| 32 | + energy = (energy - np.min(energy)) / (kJ / mol) |
| 33 | + print('E_vinyl - E_aldehyde = {:7.3f} kJ/mol'.format(energy[1] - energy[0])) |
| 34 | + |
| 35 | + # generate initial structures using metadynamics |
| 36 | + plumed_str = PLUMED_INPUT |
| 37 | + plumed_str += 'METAD ARG=CV PACE=5 SIGMA=0.25 HEIGHT=5\n' |
| 38 | + metadynamics = Metadynamics(plumed_str) |
| 39 | + |
| 40 | + # create 40 identical walkers |
| 41 | + walkers = Walker( |
| 42 | + aldehyd, |
| 43 | + hamiltonian=mace, |
| 44 | + temperature=300, |
| 45 | + metadynamics=metadynamics, |
| 46 | + ).multiply(4) |
| 47 | + |
| 48 | + # do MTD and create large dataset from all trajectories |
| 49 | + outputs = sample(walkers, steps=2000, step=20, start=1000) |
| 50 | + data_mtd = sum([o.trajectory for o in outputs], start=Dataset([])) |
| 51 | + |
| 52 | + # initialize walkers for umbrella sampling |
| 53 | + walkers = [] |
| 54 | + for i, center in enumerate(np.linspace(1, 3, num=16)): |
| 55 | + bias = get_bias(kappa=1500, center=center) |
| 56 | + hamiltonian = mace + bias |
| 57 | + walker = Walker(alcohol, hamiltonian=hamiltonian, temperature=300) |
| 58 | + walkers.append(walker) |
| 59 | + quench(walkers, data_mtd) # make sure initial structure is reasonable |
| 60 | + replica_exchange(walkers, trial_frequency=100) # use REX for improved sampling |
| 61 | + |
| 62 | + outputs = sample(walkers, steps=1000, step=10) |
| 63 | + |
| 64 | + |
| 65 | +if __name__ == '__main__': |
| 66 | + with psiflow.load() as f: |
| 67 | + main() |
0 commit comments