

Python interface
================

For different reasons, it can be useful to have directly access to the variables and subroutines in MW from outside the code itself.
To this end, we interfaced the Fortran source code with python. This means that one can run a MW simulation, exploiting the velocity of a Fortran code, from a python script which is easy to write and fast to modify.

Building the interface
----------------------
Building the interface means creating a library which contains links to most of the source code of MW. The compilation is included in the *Makefile* and is done by the command

```
make python
```

The tool used to compile the library is **F2Py**, usually given as part of the python distribution, or of the numpy package. The path to F2Py should be given in the *config.mk* file (see the [installation instructions](install)).
However, F2Py does not support the interfacing of Fortran derived types, which are widely used in MW.
To gain access to these derived types and their variables, we prepare wrappers using the software **f90wrap**, available open source [here](https://github.com/jameskermode/f90wrap), through the command

```
f90wrap -m metalwalls -k .f2py_f2cmap $(src_python)
```

The `$(src_python)` are the preprocessed source files, stored in *./build/preprocess*.
This command builds a wrapper for each Fortran file named *f90wrap_fortran_name.f90* that handles derived types and it builds an additional library *metalwalls.py* that imports the F2Py library. However, before compiling this interface, some modifications have to be made to the *f90wrap_\*.f90* files because there are some subroutines with the same name in different modules. This is done by running the script

```
./scripts/modify_python_f90wrap.sh
```

which renames the problematic subroutines by prepending the module's name to them. When adding new files, it could be necessary to change this script to handle any new conflicting subroutine. 

The Fortran-python interface is then compiled using

```
f2py -c --f90exec=mpif90 --fcompiler=gnu95 -m _metalwalls ./build/src/*.o --f90flags="$(F90FLAGS)" f90wrap_*.f90
```

which creates a library with a name resembling *_metalwalls.cpython-36m-x86+64-linux-gnu.so*.

All files generated by the compilation are placed in the folder *./build/python* and should be included in the $PYTHONPATH so that the libraries can be imported. This is easily done by typing

``` 
export PYTHONPATH=$PYTHONPATH:path_to_python_sources
```

or adding this line to the *~/.bashrc* file.

The interface contains wrapper functions to all module variables and subroutines present in the MW source files, except for a few exceptions. Indeed, f90wrap allows to handle derived types but not arrays of derived types. In MW, there are only a few such arrays, naming the arrays `ions(:)`, `electrodes(:)`, `molecules(:)` and `thomas_fermi(:)` in the `MW_system` module, which cannot therefore be accessed. If needed, new Fortran subroutines have to be especially added. For these and other reasons, the compilation takes into account all source files except *electrode.f90*, *ion.f90*, *molecule.f90*, *thomas_fermi.f90*,*localwork.f90*, *kinds.F90* and *constants.F90*.

Running MW through the python interface
-------------------
Once imported, the library gives access to both types and subroutines defined in the Fortran modules. To use a variable, it first has to be defined, for example to define a variable `my_system` of type `MW_system_t` defined in the module `MW_system`, the command is
```python
my_system = metalwalls.mw_system.MW_system_t()
```
Similarly, to use the subroutine that runs the simulation steps defined in the module `MW_tools` the command is
```python
metalwalls.mw_tools.run_step(my_system, my_parallel, my_algorithms, istep, num_steps, do_output, step_output_frequency, True, True)
```

To obtain information on a variable or how to use a specific subroutine, type in the prompt
```python
>>> print(my_system.__doc__)

        Type(name=mw_system_t)
        
        
        Defined at ./src/system.F90 lines 114-262
        
        
>>> print(metalwalls.mw_tools.__doc__)

    Module mw_tools
    
    
    Defined at ./src/tools.F90 lines 1-761
    
    
>>> print(metalwalls.mw_tools.run_step.__doc__)

        run_step(self, parallel, algorithms, istep, num_steps, do_output, step_output_frequency, pre, post)
        
        
        Defined at ./src/tools.F90 lines 395-433
        
        Parameters
        ----------
        system : Mw_System_T
        parallel : Mw_Parallel_T
        algorithms : Mw_Algorithms_T
        istep : int
        num_steps : int
        do_output : bool
        step_output_frequency : int
        pre : bool
        post : bool
    

>>> print(metalwalls.mw_tools.initialize.__doc__)

        do_output = initialize(self, parallel, algorithms)
        
        
        Defined at ./src/tools.F90 lines 31-114
        
        Parameters
        ----------
        system : Mw_System_T
        parallel : Mw_Parallel_T
        algorithms : Mw_Algorithms_T
        
        Returns
        -------
        do_output : bool
        

```

Using MW through the python interface in general requires a deep knowledge of how the code works and is built. An additional library with several utilities was written to allow for simple operations. These python scripts are in the *./src/* folder and copied to the *./build/python* folder upon compilation:

- *simulation.py* defines basic python functions to initialize, setup, run and finalize a simulation
- *piston.py* allows to run simulations with electrodes used as pistons, when electrodes are formed of several smaller electrodes (with different species)
- *thermodynamic_integration.py* defines python functions to run thermodynamic integrations between initial and final states with different charges, Lennard-Jones parameters or introducing a soft-core Lennard-Jones particle. 
- *mw.py* imports all other scripts

These scripts are made for specific simulations and should be modified and adapted with care.


Step-by-step script explanation
------------------------------

To understand how the python interface can be used, we will go through a script that computes a free energy difference using thermodynamic integration for an aqueous system of 299 water molecules and one ION. This example is given in the folder *./example/python-thermo_integration_lj_lj_brute_force*. In the initial state, the Lennard-Jones parameters of the ION are those of a sodium ion whereas in the final state they are those of a potassium ion. Therefore we computed the difference in hydration free energy between the sodium and potassium ions.

By naming the Hamiltonians corresponding to the initial and final state respectively $`H_{Na}`$ and $`H_{K}`$, and define an hybrid Hamiltonian $`H(\lambda) = (1-\lambda) H_{Na} + \lambda H_K`$, the thermodynamic integration states that

```math
\Delta F = \int \limits _0 ^1 \mathrm{d} \lambda \left< \frac{\partial H}{\partial \lambda} \right>_\lambda
```

This thus requires the computation of the average of $`\frac{\partial H}{\partial \lambda} = H_K - H_{Na}`$ along a simulation run with the hybrid Hamiltonian $`H(\lambda)`$. This kind of simulation cannot be directly run with MW and would need significant modifications in the source code and a complex input file. Here is how to do this thermodynamic integration using the python interface, modifying the forces and computing the energy derivative using MW subroutines.

The first element in the python script must be to import the MW library, that also needs the **numpy** package

```python
import numpy as np
import metalwalls
```

Next, we need to initialize the MPI parallel environment. To this end, the python distribution should have the **mpi4py** package installed. First import the library

```python
me = 0

from mpi4py import MPI
comm = MPI.COMM_WORLD         
fcomm = comm.py2f()           
me = comm.Get_rank()        # rank of current processor
nprocs = comm.Get_size()    # total number of processors 
```

Then create a variable `my_parallel` as the derived type `MW_parallel_t` defined in the module `mw_parallel`

```python
my_parallel = metalwalls.mw_parallel.MW_parallel_t()
```

Start the parallel environment using the subroutine `start_parallel_python` defined in the module `mw_tools`. The communicator to be given to this subroutine is not directly the one obtained from mpi4py, but is obtained doing `fcomm = comm.py2f()`

```python
metalwalls.mw_tools.start_parallel_python(my_parallel, fcomm)
```

Next, loop over different values of `lambd` to compute the energy derivative needed for the thermodynamic integration

```python
lambdas = np.arange(0, 1.1, 0.1)         # list of lambda values
energy_derivative = np.zeros((11, 3))    # array 0:lambda, 1:dU, 2:std error
for i, lambd in enumerate(lambdas):
   energy_derivative[i, 0] = lambd
   energy_derivative[i, 1], energy_derivative[i, 2] = run_lambda(lambd, 100, 50)
   if me == 0: print('Lambda {:4.2f}: dU {} hartree'.format(lambd, energy_derivative[i, 1])
   comm.Barrier()
```

The free energy difference is then obtained by integrating the energy derivative

```python
ti = np.trapz(energy_derivative[:, 1], x=energy_derivative[:, 0])  # in hartree
```

Last, finalize the parallel environment using the subroutine `finalize_parallel` 

```python
metalwalls.mw_tools.finalize_parallel()
```

For each `lambd` value, a simulation with a hybrid Hamiltonian is run, using the python function `run_lambda` that manipulates the forces and extracts the energies

```python
def run_lambda(lambd, num_steps_per_cycle, num_cycles):
```

The steps to run a simulation follow the same steps as in the *main.f90*. First define essential variables

```python
   my_system = metalwalls.mw_system.MW_system_t()
   my_algorithms = metalwalls.mw_algorithms.MW_algorithms_t()
```

The `system` and `algorithms` variables in MW, along with the `parallel` variable are the three most important variables (all derived types) in MW that contain all the necessary information about the simulation.

Then initialize the simulation

```python
   #Initialize simulation, read input files
   do_output = metalwalls.mw_tools.initialize(my_system, my_parallel, my_algorithms)
```

This step reads the input files (*runtime.inpt* and *data.inpt*), allocates arrays, sets up various flags and quantities, distributes the work amongst processes and outputs the simulation informations in the *run.out*.

Since the thermodynamic integration is between two sets of Lennard-Jones parameters (one for $`Na^+`$ and one for $`K^+`$), we should define those two sets. To this end, first copy the variable `my_system.lj` created when reading the input files, then modify the desired parameters

```python
   # Create two lj derived types for Na and K
   ljNa = metalwalls.mw_lennard_jones.copy_type(my_system.lj)
   ljK = metalwalls.mw_lennard_jones.copy_type(my_system.lj)
    
   # ION-O parameters
   epsNa = 0.521578    # kJ/mol
   sigmaNa = 2.876500  # angstroms
   epsK = 0.5216       # kJ/mol
   sigmaK = 3.25       # angstroms
   
   # Modify the ION-O parameters: O=1, ION=4 (Fortran indexes)
   # (no ION-ION parameter since there is only one ION)
   metalwalls.mw_lennard_jones.set_parameters(ljNa, 1, 4, epsNa, sigmaNa)
   metalwalls.mw_lennard_jones.set_parameters(ljK, 1, 4, epsK, sigmaK)
```

Then setup the simulation for step 0, that is, precompute electrode charges, dipoles, forces

```python
   #Setup simulation
   step_output_frequency = metalwalls.mw_tools.setup(my_system, my_parallel, my_algorithms, do_output)
```

And initialize various timers

```python
   #Initialize run, timers 
   metalwalls.mw_tools.run_step(my_system, my_parallel, my_algorithms, 0, 0, do_output, step_output_frequency, True, False)
```

The 4th parameter is the number of steps (here 0) and the last parameters are respectively whether to start or stop the timers. 

Next, loop over the number of steps of the simulation. In this case, there are `num_cycles` cycles of `num_steps_per_cycle` steps each. This way, the energy derivative is computed only once per cycle.

```python
   # Loop over the cycles
   energy_diff = np.zeros(num_cycles)
   for cycle in range(num_cycles):
      ## Compute energy derivative dU/dlambd = U2 - U1
      # Compute the energy with the Na parameters
        # set all energies to 0
      metalwalls.mw_lookup_table.lut_set_all(my_system.energy, 0.0)
        # set the system.lj parameter to the Na parameters
      my_system.lj = ljNa
        # compute the Lennard-Jones energy
      metalwalls.mw_system.energy_lj(my_system)
        # get the Lennard-Jones energies from the lookup table
      hNa = metalwalls.mw_lookup_table.lut_get(my_system.energy, "Lennard-Jones")
      hNa += metalwalls.mw_lookup_table.lut_get(my_system.energy, "Lennard-Jones_cutoff")
        
      # Compute the energy with the K parameters
        # set all energies to 0
      metalwalls.mw_lookup_table.lut_set_all(my_system.energy, 0.0)
        # set the system.lj parameter to the K parameters
      my_system.lj = ljK
        # compute the Lennard-Jones energy
      metalwalls.mw_system.energy_lj(my_system)
        # get the Lennard-Jones energies from the lookup table
      hK = metalwalls.mw_lookup_table.lut_get(my_system.energy, "Lennard-Jones")
      hK += metalwalls.mw_lookup_table.lut_get(my_system.energy, "Lennard-Jones_cutoff")
        
      # Compute the energy difference
      energy_diff[cycle] = hK - hNa
```

Since the hybrid Hamiltonian is a linear combination of the Na and K Hamiltonians, the energy derivative is the difference between the energies. The energies in MW are stored in a lookup table, set up in the *lookup_table.f90* file, and different subroutines to set or get energies have to be used to modify them.

Then the simulation has to be run for `num_steps_per_cycle` using hybrid forces. A step is made of 3 sub-steps: first a setup (where the positions of the atoms are updated and the electrode charges, dipoles or radius are computed if needed), then the new forces are computed and lastly an update step (where the velocities are updated - second step of the Velocity Verlet - and information is output). In this case, the computation of forces needs to be modified

```python
      # Run steps
      for step in range(num_steps_per_cycle):
         ## First step: setup the step, first part of time integration
         metalwalls.mw_tools.step_setup(my_system, my_algorithms, my_parallel, (cycle + 1) * num_steps_per_cycle + step)
            
         ## Second step: compute forces
         ## Usually done by 
         ## metalwalls.mw_tools.step_compute_forces(my_system, my_parallel, (cycle + 1) * num_steps_per_cycle + step)
           # Compute intramolecular forces
         my_system.intramolecular_forces_ions[:,:,:] = 0.0
           # Compute coulomb forces
         my_system.coulomb_forces_ions[:,:,:] = 0.0
         metalwalls.mw_system.forces_compute_coulomb(my_system)
           # Compute intermolecular forces
         uNa = np.zeros(np.shape(my_system.intermolecular_forces_ions))
         uK = np.zeros(np.shape(my_system.intermolecular_forces_ions))
             # compute with model Na and store in uNa
         my_system.intermolecular_forces_ions[:,:,:] = 0.0
         my_system.lj = ljNa
         metalwalls.mw_system.forces_compute_lj(my_system)
         uNa[:, :, :] = my_system.intermolecular_forces_ions[:, :, :]
             # compute with model K and store in uK
         my_system.intermolecular_forces_ions[:,:,:] = 0.0
         my_system.lj = ljK
         metalwalls.mw_system.forces_compute_lj(my_system)
         uK[:, :, :] = my_system.intermolecular_forces_ions[:, :, :]
             # the forces should be the linear combination of uNa and uK
         my_system.intermolecular_forces_ions[:, :, :] = (1 - lambd) * uNa[:, :, :] + lambd * uK[:, :, :]
           # Gather all forces (coulomb, intermolecular, intramolecular)
         metalwalls.mw_system.forces_gather(my_system)
        
         ## Third step: second part of time integration, outputs
         metalwalls.mw_tools.step_update_output(my_system, my_algorithms, (cycle + 1) * num_steps_per_cycle + step, do_output, step_output_frequency)
```

After all steps have been computed, stop the timers and then finalize the simulation without finalizing the parallel environment (it can be initialized only once per python execution)

```python
   #Finalize run, timers
   metalwalls.mw_tools.run_step(my_system, my_parallel, my_algorithms, num_steps_per_cycle * num_cycles, 0, do_output, step_output_frequency, False, True)

   #Finalize   
   metalwalls.mw_tools.finalize(my_system, my_parallel, my_algorithms)

   mean_dU = np.mean(energy_diff)
   err_dU = np.std(energy_diff) / num_cycles**0.5

   return mean_dU, err_dU
```

A more straightforward way is implemented in the *thermodynamic_integration.py* script, used in the folder *./example/python-thermo_integration_lj_lj*, that does not compute the linear combination of forces but uses a combination of the parameters and a single force computation to obtain the same result. 