Build 2+3-body Mapped GP

After the on-the-fly training is complete, we can play with the force field we obtained. We are going to do the following things:

1. Parse the on-the-fly training trajectory to collect training data

2. Reconstruct the GP model from the training trajectory

3. Build up Mapped GP (MGP) for accelerated force field, and save coefficient file for LAMMPS

4. Use LAMMPS to run fast simulation using MGP pair style

Parse OTF log file

After the on-the-fly training is complete, we have a log file and can use the otf_parser module to parse the trajectory.

import numpy as np
from flare import otf_parser

logdir = '../../../tests/test_files'
file_name = f'{logdir}/AgI_snippet.out'
hyp_no = 2 # use the hyperparameters from the 2nd training step
otf_object = otf_parser.OtfAnalysis(file_name)

Construct GP model from log file

We can reconstruct GP model from the parsed log file (the on-the-fly training trajectory). Here we build up the GP model with 2+3 body kernel from the on-the-fly log file.

gp_model = otf_object.make_gp(hyp_no=hyp_no)
gp_model.parallel = True
gp_model.hyp_labels = ['sig2', 'ls2', 'sig3', 'ls3', 'noise']

# write model to a binary file
gp_model.write_model('AgI.gp', format='json')

The last step write_model is to write this GP model into a binary file, so next time we can directly load the model from the pickle file as

from flare.gp import GaussianProcess

gp_model = GaussianProcess.from_file('AgI.gp.json')

Map the GP force field & Dump LAMMPS coefficient file

To use the trained force field with accelerated version MGP, or in LAMMPS, we need to build MGP from GP model. Since 2-body and 3-body are both included, we need to set up the number of grid points for 2-body and 3-body in grid_params. We build up energy mapping, thus set map_force=False. See MGP tutorial for more explanation of the MGP settings.

from flare.mgp import MappedGaussianProcess

grid_params = {'twobody':   {'grid_num': [64]},
               'threebody': {'grid_num': [20, 20, 20]}}

data = gp_model.training_statistics
lammps_location = 'AgI_Molten'

mgp_model = MappedGaussianProcess(grid_params, data['species'],
    var_map=None, lmp_file_name='AgI_Molten', n_cpus=1)
mgp_model.build_map(gp_model)

/Users/yuxie/opt/anaconda3/lib/python3.8/site-packages/flare/mgp/mapxb.py:519: UserWarning: The minimal distance in training data is lower than the current lower bound, will reset lower bound to 2.129780094032889
  warnings.warn(
/Users/yuxie/opt/anaconda3/lib/python3.8/site-packages/flare/mgp/mapxb.py:519: UserWarning: The minimal distance in training data is lower than the current lower bound, will reset lower bound to 2.129780094032889
  warnings.warn(
/Users/yuxie/opt/anaconda3/lib/python3.8/site-packages/flare/mgp/mapxb.py:519: UserWarning: The minimal distance in training data is lower than the current lower bound, will reset lower bound to 2.129780094032889
  warnings.warn(

The coefficient file for LAMMPS mgp pair_style is automatically saved once the mapping is done. Saved as lmp_file_name.

Run LAMMPS with MGP pair style

With the above coefficient file, we can run LAMMPS simulation with the mgp pair style. First download our mgp pair style files, compile your lammps executable with mgp pair style following our instruction in the Installation section.

1. One way to use it is running lmp_executable < in.lammps > log.lammps with the executable provided in our repository. When creating the input file, please note to set

newton off
pair_style mgp
pair_coeff * * <lmp_file_name> <chemical_symbols> yes/no yes/no

An example is using coefficient file AgI_Molten.mgp for AgI system, with two-body (the 1st yes) together with three-body (the 2nd yes).

pair_coeff * * AgI_Molten.mgp Ag I yes yes

2. Another way is to use the ASE LAMMPS interface

import os
from flare.utils.element_coder import _Z_to_mass, _element_to_Z
from flare.ase.calculator import FLARE_Calculator
from ase.calculators.lammpsrun import LAMMPS

from ase import Atoms

# create test structure
species = otf_object.gp_species_list[-1]
positions = otf_object.position_list[-1]
forces = otf_object.force_list[-1]
otf_cell = otf_object.header['cell']
structure = Atoms(symbols=species, cell=otf_cell, positions=positions)

# get chemical symbols, masses etc.
species = gp_model.training_statistics['species']
specie_symbol_list = " ".join(species)
masses=[f"{i} {_Z_to_mass[_element_to_Z[species[i]]]}" for i in range(len(species))]

# set up input params
parameters = {'command': os.environ.get('lmp'), # set up executable for ASE
              'newton': 'off',
              'pair_style': 'mgp',
              'pair_coeff': [f'* * {lammps_location + ".mgp"} {specie_symbol_list} yes yes'],
              'mass': masses}
files = [lammps_location + ".mgp"]

# create ASE calc
lmp_calc = LAMMPS(label=f'tmp_AgI', keep_tmp_files=True, tmp_dir='./tmp/',
        parameters=parameters, files=files, specorder=species)

structure.calc = lmp_calc

# To compute energy, forces and stress
# energy = structure.get_potential_energy()
# forces = structure.get_forces()
# stress = structure.get_stress()

/Users/yuxie/opt/anaconda3/lib/python3.8/site-packages/ase/calculators/lammpsrun.py:191: UserWarning: You are using an old syntax to set 'parameters'.
Please use LAMMPSRUN.set().
  warnings.warn(self.legacy_warn_string.format("parameters"))

3. The third way to run LAMMPS is using our LAMMPS interface, please set the environment variable $lmp to the executable.

from flare import struc
from flare.lammps import lammps_calculator

# lmp coef file is automatically written now every time MGP is constructed

# create test structure
species = otf_object.gp_species_list[-1]
positions = otf_object.position_list[-1]
forces = otf_object.force_list[-1]
otf_cell = otf_object.header['cell']
structure = struc.Structure(otf_cell, species, positions)

atom_types = [1, 2]
atom_masses = [108, 127]
atom_species = [1, 2] * 27

# create data file
data_file_name = 'tmp.data'
data_text = lammps_calculator.lammps_dat(structure, atom_types,
                                         atom_masses, atom_species)
lammps_calculator.write_text(data_file_name, data_text)

# create lammps input
style_string = 'mgp'
coeff_string = '* * {} Ag I yes yes'.format(lammps_location)
lammps_executable = '$lmp'
dump_file_name = 'tmp.dump'
input_file_name = 'tmp.in'
output_file_name = 'tmp.out'
input_text = \
    lammps_calculator.generic_lammps_input(data_file_name, style_string,
                                           coeff_string, dump_file_name)
lammps_calculator.write_text(input_file_name, input_text)

# To run lammps and get forces
# lammps_calculator.run_lammps(lammps_executable, input_file_name,
#                              output_file_name)

# lammps_forces = lammps_calculator.lammps_parser(dump_file_name)