Tutorial: Add a Model System

It is very easy and straightforward to add your own model system to the PySurf Framework via the Plugin engine. In the models, specific user input can be put in and your model is natively added to the program package. Just put your python script in the plugin folder.

Example: 1D Harmonic Oscillator

In this tutorial, we will implement a simple model with one or two harmonic oscillators.The user can decide, whether the second surface is included and whether it should be shifted, as well as the frequencies and the offset energy.

Step 0: Colt Factory

Models are written as Plugins of the Model Colt-Factory, i.e. they are inherited from the Factory.

1from pysurf.spp import Model
2from pysurf.system import Mode

Step 1: Basic Features

Every Model has to provide some features:

implemented

which properties are provided by the model, e.g. [‘energy’, ‘gradient’, ‘fosc’]

modes

list of the ground state modes for the Samplers

@classmethod def from_config

takes user input if questions are there and has to call the init function

def get(self, request):

fill the request with the demanded data

 1 class HarmonicOscillator1D():
 2     """ Model for a 1D harmonic oscillator with 1 or 2 PES """
 3
 4     implemented = ['energy', 'gradient']
 5     crd = [0.0]
 6     frequencies = [1.0]
 7     displacements = [[1.0]]
 8     modes = [Mode(freq, dis) for freq, dis in zip(frequencies, displacements)]
 9
10     @classmethod
11     def from_config(cls, config):
12         """initialize class with a given questions config!"""
13
14     def get(self, request):
15         """get requested info"""

Step 2: User Input

PySurf is build around the Colt framework, developed along the lines of this project. To specify specific input needed for your class you simply use the _questions string: In our example we will need user input for the shape of the potentials:

  • E0/E1

    offset energy of the corresponding potential

  • x0/x1

    shift of the harmonic oscillator in x direction

  • w0/w1

    frequency of the corresponding PES

  • npes

    number of PES. The parameters for the second PES are only asked if npes = 2

 1 class HarmonicOscillator1D(Model):
 2 """ Model for a 1D harmonic oscillator with 1 or 2 potential energy surfaces """
 3
 4 _questions = """
 5 e0 = 0.0 :: float
 6 w0 = 1.0 :: float
 7 x0 = 0.0 :: float
 8 # Number of potential energy surfaces
 9 npes = 1 :: str ::
10
11 [npes(1)]
12
13 [npes(2)]
14 e1 = 1.0 :: float
15 w1 = 1.0 :: float
16 x1 = 1.0 :: float
17 """
18
19 implemented = ["energy", "gradient"]
20 masses = [1.0]
21 crd = [0.0]
22 frequencies = [1.0]
23 displacements = [[1.0]]
24 modes = [Mode(freq, dis) for freq, dis in zip(frequencies, displacements)]
25
26 @classmethod
27 def from_config(cls, config):
28     e0 = config['e0']
29     w0 = config['w0']
30     x0 = config['x0']
31     npes = int(config['npes'].value)
32
33     config_npes = config['npes']
34     return cls(e0, w0, x0, npes, config_npes)
35
36 def __init__(self, e0, w0, x0, npes, config_npes):
37     self.frequencies = [w0]
38     self.crd = [x0]
39     self.npes = int(npes)
40     self.w = [w0]
41     self.x = [x0]
42     self.e = [e0]
43
44     if self.npes == 2:
45         self.w += [config_npes['w1']]
46         self.x += [config_npes['x1']]
47         self.e += [config_npes['e1']]

According to the Colt style, all arguments in the main question block are given in the init method explicitely, whereas other question blocks are passed as config of the block.

Step 3: Implement _energy and _gradient function

The next step is to implement the actual model, i.e. the properties of the model. In our case we want to provide the energies and the gradients of the surfaces and thus private methods for _energy and _gradient are implemented.

 1 class HarmonicOscillator1D(Model):
 2 """ Model for a 1D harmonic oscillator with 1 or 2 potential energy surfaces """
 3
 4   ...
 5
 6     def _energy(self, x):
 7         energy = []
 8         for i in range(self.npes):
 9             energy += [0.5*self.w[i]*(x - self.x[i])**2 + self.e[i]]
10         energy = np.array(energy).flatten()
11         return energy
12
13     def _gradient(self, x):
14         gradient = {}
15         for i in range(self.npes):
16             gradient[i] = np.array(self.w[i]*(x - self.x[i]))
17         return gradient

The _energy function takes a coordinate position, i.e. a numpy array and returns an array with one or two entries, depending on whether the surface contains one or two potential energy surfaces. The energy is calculated according to the formula $0.5*omega*(x-x0)^2$

The _gradient function takes a coordinate position, i.e. a numpy array and returns a dictionary. The keys in the dictionary are the state numbers as integers, i.e. 1 or 2 and the values are the gradients of the corresponding state. In our specific case, such a dictionary may look like: {1: [0.5], 2: [0.2]}

Step 4: Implement the get function

The get function is called with the request as parameter. It has to fill in the desired results from the model into the request instance.

 1 class HarmonicOscillator1D(Model):
 2 """ Model for a 1D harmonic oscillator with 1 or 2 potential energy surfaces """
 3
 4 ...
 5
 6 def get(self, request):
 7     """the get function returns the adiabatic energies as well as the
 8        gradient at the given position crd. Additionally the masses
 9        of the normal modes are returned for the kinetic Hamiltonian.
10     """
11     crd = request.crd
12     print('crd', crd)
13     for prop in request:
14         if prop == 'energy':
15             request.set('energy', self._energy(crd))
16         if prop == 'gradient':
17             request.set('gradient', self._gradient(crd))
18     return request

The get function checks which properties are demanded in the request and fills them into the request, using the _energy and _gradient functions.