Procedures

Insert a VdW screening pair in the link atom structure

Create adjacency matrix $\mathbb C_1$ from connectivity lists.
Array i12 and n12 contain the connectivity list in the following format: i12(0:n(j),j) contains the index of all the atoms connected to atom with index j.

Compute the torsional angle between four atoms specified with indices idx; results are in range [-pi;pi]

Trigonal center

Initialize arrays used in calculation of angle bending functions

Compute angle-bending terms of the potential energy function.
Simple angle terms are computed according to the formula: $$U_{angle} = \sum_i k_i \Delta \theta_i^2 \large(1 + \sum_{j=1}^4 k^{(j+2)} \Delta \theta_i^j \large)$$ $$\Delta \theta_i = \theta_i - \theta^{(eq)}_i$$
Out-of plane angle are more complex. First, central atom has to be a trigonal center, the other two atoms together with the auxliary atom (that is the remaining one connected to the trigonal center) define the projection plane. During the first run the auxiliary atom is found and saved. Then, the trigonal center is projected on the plane defined by the other three atoms, and the angle is the one defined by the projection (which is the vertex, and the other two atoms -- the auxiliary is excluded). Then the same formula used for simple angle terms is used. Trigonal center Normal vector of the projection plane

Initialize angle-torsion coupling potential arrays

If there are no bonds in the system just return, there is nothing to do.

No parameters are defined, nothing to do.

Initialize array used in calculation of bond stratching terms of potential energy

Compute the bond-stretching terms of the potential energy.
They are computed according to the general formula adopted in AMOEBA Force Field: $$U_{bond} = \sum_i k_i \Delta l_i^2 \large(1 + k^{(3)}\Delta l_i + k^{(4)}\Delta l_i^2 \large)$$ $$\Delta l_i = l_i - l^{(eq)}_i$$

Just terminate every "submodule" in bonded, deallocating arrays and disabling the potential terms

Build connectivity matrices up to $\mathbb C_n$ starting from $\mathbb C_1$. Results are stored in an array of boolean sparse matrix in Yale format in such a way that $res(n) := \mathbb C_n$; since FORTRAN is 1-based the useless $\mathbb C_0$ is not stored.

Builds the adjacency matrix of polarization groups starting from atomic adjacency matrix and list of polarization groups indices.

Convert a string coming from C into a Fortran string

Return true if the current forcefield is AMOEBA, and false in all other cases.

Return the c-pointer to the array containing the coordinates of MM atoms.

Return the c-pointer to the array containing the coordinates of MM atoms.

Return the c-pointer to the array containing the coordinates of polarizable atoms.

Return the c-pointer to the array containing the induced dipoles on polarizable sites.

Return the number of MM atoms in the system.

Return the number of dipole's set for the current Force-Field.

Return the number of polarizable atoms in the system.

Return the c-pointer to the array containing the map from polarizable to MM atoms.

Return the c-pointer to the array containing the static source of the electrostatic field.

Return the c-pointer to the array containing the coordinates of MM atoms.

Initalize OMMP System Object from .mmp file

Initialize the library using a Tinker xyz and a Tinker prm

Print a summary of the system input on standard output.

Print a summary of the system input on file.

use mod_mmpol, only: pol_atoms

Set the verbosity level of the library to verb

Terminate a OMMP System Object

Check if adjacency matrix up to nth order is present in topology object. If it is not present, update the topology accordingly.

Check if n new screening pairs could be allocated in la structure. If the allocated arrays are too small, they are reallocated on-the-fly.

Handles the memory errors (including soft limit) during memory allocation

Handles the memory errors (including soft limits) during the deallocation

This subroutine changes the output file for mmpol to a file defined by filename. A file has already been set, close it before proceed.

Compress the data in uc_data to the same Yale sparse format described in s

This subroutine takes as input a sparse matrix (rank [n]) in an uncompressed yale format [uc_list], as a rectangular matrix ([n] x max_el_per_row) and the actual number of items [nit] for each row (remaining elements are not considered) and compress in a Yale format sparse matrix [s]. The task is parallelized to handle large matrices.

Evaluate the z value at position (x, y) of a surface built as a bicubic spline interpolating the points
(xgrd($x_i$, $y_i$), ygrd($x_i$, $y_i$), vxgrd($x_i$, $y_i$)).
In order to do so, also the derivatives of the surface at the points of the grid are needed; in particular if we consider the surface points as $$(x_i, y_i, V(x_i, y_i))$$ also value of $\frac{\partial V}{\partial x}$, $\frac{\partial V}{\partial y}$ and $\frac{\partial^2 V}{\partial x \partial y}$ at grid points are needed.

Conjugate gradient solver (TODO)

Compute the cosine of torsional angle between four atoms specified with indices idx

This function compute the coulomb kernel for the distance vector dr and its derivative up to the value required by maxder.

Count the number of occurence of substring c in string s, and return the number of occurence, if c is not contained in s, zero is returned.

Create a new bond between atoms i and j in the topology.

Explicitly construct polarization tensor in memory. This routine is only used to accumulate results from dipole_T and shape it in the correct way.

Computes the cross product between two vectors of dimension 3. Used as an inlinable function in geometric manipulations

Compute the cyclic interpolating cubic spline (2D) that passes for
points $(x_i, y_i)$. Each segment is described by the curve: $$S_i(x) = a_i + b_i(x-x_i) + c_i (x-x_i)^2 + d_i(x-x_i)^3$$ The algorithm used to compute the coefficients is taken from "Numerical Algorithm with C" - Gisela ENGELN-MULLGES Frank UHLIG (10.1007/978-3-642-61074-5)

Output a 1D-matrix of real in a well formatted way

Output a 2D-matrix of real in a well formatted way

This subroutine computes the damped coulomb kernel between two atoms. Note that this only makes sense between two MM atoms, as it is only used to compute the field that induces the point dipoles!

perform Pulay's direct inversion in the iterative subspace extrapolation:

This subroutine compute the interaction tensor (rank 3) between two polarizable sites i and j. This tensor is built according to the following rules: ... TODO

Computes the electric field of a trial set of induced point dipoles at polarizable sites. This is intended to be used as matrix-vector routine in the solution of the linear system.

Computes the electric field of static multipoles at induced dipoles sites. This is only intended to be used to build the RHS of the linear system. This field is modified by the indroduction of the damped kernels and by the connectivity-based screening rules.

Computes the electric potential, field and field gradients of static multipoles at all sites (polarizable sites are a subset of static ones)

Computes the electrostatic quantities (that means nuclear-MM
interaction terms) needed to perform an energy calculation. Computes: (1) EF of nuclei at polarizable sites (2) V of the whole MM system at QM sites

Computes the electrostatic quantities (that means nuclear-MM
interaction terms) needed to perform a gradient calculation. Computes: (1) GEF on nuclei at polarizable sites (2) EF, GEF, HEF of nuclei at static sites (3) EF of whole MM system at QM sites

This function computes the interaction energy of static electric multipoles

This function computes the interaction energy of static electric multipoles

Prints a message and exit from the program. This function should be used in all the conditions where the program cannot proceed.

This subroutine computes the potential generated by the static multipoles to a set of arbitrary coordinates, without applying any screening rules.

Computes the electric field of a trial set of induced point dipoles at polarizable sites. This is intended to be used as matrix-vector routine in the solution of the linear system.

This subroutine computes the potential generated by the static multipoles to a set of arbitrary coordinates, without applying any screening rules.

Once that the neighbor list have been initialized and updated, this function provide a logical array for atom i with all interactions that should be computed and corresponding distances.

This function is intended to check if the ff described by prm_type is AMOEBA (or amoeba-like) or AMBER or FF of another kind. A FF is considered to be AMOEBA if: it contains multipole keywords (and no charge keywords) and polarization MUTUAL. A FF is considered do be AMBER if contains charge keyword and no multipole keyword.

This subroutine guess the connectivity of the system from the coordinates and the atomic number of its atoms. It is based on distances and atomic radii, so it can easily fail on distorted geometries. It should be used only when the the bonds of the molecule are not availble in any other way; it is often used to assign connectivity to a QM part that does not have any. H B C N O F Si P S Cl As Se Br Te I

Allocate a 1-dimensional array of integers

Allocate a 2-dimensional array of integers

Allocate a 3-dimensional array of integers

Free a 1-dimensional array of integers

Free a 2-dimensional array of integers

Free a 3-dimensional array of integers

Compute torsion potential

Initialize improper torsion potential arrays

Compute torsion potential

Insert in the bonded parameter required for the link atom between iqm and imm Link atoms should never appear in bonded interactions!

Initialize the quantities needed for vdw screening due to the presence of a link atom between iqm and imm

Computes the Jacobian matrix for the inplane angle definition. It computes the Jacobian for the normal angle using the projected point (R) as central point. Then projects $J_r = \frac{\partial \theta}{\partial \vec{R}}$ onto A, B, C, and X (auxiliary point). The projection is done computing the 3x3 matrices of partial derivative of $\vec{R}$ wrt any actual point and using them to project $J_r$. [\frac{\partial \vec{R}}{\partial \vec{A}} = \begin{bmatrix} \frac{\partial \vec{R}_x}{\partial \vec{A}_x} & \frac{\partial \vec{R}_y}{\partial \vec{A}_x} & \frac{\partial \vec{R}_z}{\partial \vec{A}_x} \ !! \frac{\partial \vec{R}_x}{\partial \vec{A}_y} &

Solve the linear system directly inverting the matrix: $$\mathbf A \mathbf x = \mathbf B$$ $$\mathbf x = \mathbf A ^-1 \mathbf B$$ This is highly unefficient and should only be used for testing other methods of solution.

Decide if a string can be interpreted as an integer or not

Decide if a string can be interpreted as a real

Routine to perform matrix-vector product

Allocate a 1-dimensional array of reals

Allocate a 2-dimensional array of reals

Free a 1-dimensional array of integers

Free a 2-dimensional array of integers

This function attempt some magic to speed up the read of a large text file in order to transfer it to memory.

Compute the cartesian coordinates of link atom idx at the current geometry.

Update merged topology in linkatom object so that its coordinates are the same of mmtop and qmtop.

assemble the DIIS B matrix:

Interface to perform memory allocation within the openMMPol library, it can be called for 1,2 and 3-dimensional arrays of either integer or real

Performs the operation $res := sp1 \land \neg sp2$ on boolean sparse matrices.

Performs the operation $res := sp1 \cdot sp2$ on boolean sparse matrices; product correspond to logical and while sum correspond to logical or.
This subroutine is the one actually used in the code; It performs matrix product exploiting the sparsity, and therefore with a scaling $\mathcal O(n)$.

Performs the operation $res := sp1 \cdot sp2$ on boolean sparse matrices; product correspond to logical and while sum correspond to logical or.
This subroutine is just for test pourposes as it performs the matrix product without exploiting the sparsity, and therefore with an unfavorable scaling $\mathcal O(n^2)$.

Copy boolean sparse matrix in yale format f to t.

Deallocate a boolean sparse matrix (in Yale format)

Return the current value of maximum used memory, if an argument is present set max_used to that value, otherwise max_used is reset to -1

Routine used to initialize the memory module. It should be called during the module initialization.

Merge topologies top1 and top2 to create top3 no link between the two topologies are created. map13 and map23 are arrays mapping the atoms of top1 to top3 and the ones of top2 to top3 respectively

Interface to perform memory deallocation within the openMMPol library, it can be called for 1,2 and 3-dimensional arrays of either integer or real

Performs all the memory allocation and vector initialization needed to run the openMMPol library

Enable nonbonded part of pontential

This function read a .mmp file (revision 2 and 3) are supported and initialize all the quantities need to describe the environment within this library.

This function read a Tinker xyz and prm files and initialize all the quantities need to describe the environment within this library.

Enable link atom

Enable nonbonded part of pontential

Prints a complete summary of all the quantities stored in the MMPol module

Compute some derived quantities from the input that are used during the calculation. The upstream code have to provide cmm, q, pol, adjacency matrix and in the case of AMOEBA also multipoles rotation information, and polarization group information.
This routine
* compute connectivity lists from connected atoms
* invert polar_mm list creating mm_polar
* populate cpol list of coordinates
* compute factors for thole damping
* scales by 1/3 AMOEBA quadrupoles (?)
* Build list for polarization groups, compute groups connectivity
* performs multipoles rotation

Save the loaded system in mmpol format. Only the electrostatic part is saved, everything else is just ignored. Version 2 and 3 of the mmp format are supported, 3 is used as default.

Performs all the deallocation needed at the end of the calculation

TODO this should be improved

Solve the polarization equation for a certain external field and compute the interaction energy of the induced dipoles with themselves and fixed multipoles.

Output a message according to the verbosity level.

This function is an interface for saving an HDF5 file with all the data contained in mmpol module using mod_io

This function get an external field and solve the polarization system in the presence of the provided external field.

This is just the same as ommp_set_external_field but implicitly assuming ommp_set_external_field as false, mainly here for interface consistency with C

External interface for smartinput function

External interface for smartinput function

Takes in input a QM Helper object, with initialized atom types, and using a parameter file, it generates a OMMP System object that corresponds to the QM system. It is used for internal testing pourpose but other creative things are always possible.

Initialize arrays for out-of-plane bending potential calculation.

Computes the out-of-plane bending potential.
With Allinger formula: similarly to in plane angles, here we are considering a trigonal center, where D is the central atom and A, B, C are connected to D. Allinger formula consider the angle between vector $\vec{AD}$ and the normal vector of plane ABC, using $\frac{\pi}{2}$ as implicit equilibrium value. The formula for this potential term is: $$U_{out-of-plane} = \sum_i k_i \chi_i^2 \large(1 + \sum_{j=1}^4 k^{(j+2)} \chi_i^j \large)$$

Initialize arrays needed to compute pi-torsion potential

Compute pi-torsion terms of the potential.
This potential is defined on a $\pi$-system, and uses the coordinates of six atoms A...F the central "double" bond is A-B, then C and D are connected to A while E and F are connected to B. So two plane ACD and BEF are defined. The potential is computed using the dihedral angle of the normal vector of those two planes, connected by segment A-B ($\theta$).
The formula used is: $$U_{\pi-torsion} = \sum_i k_i \large(1 + cos(2\theta-\pi) \large)$$

Main driver for the calculation of induced dipoles. Takes electric field at induced dipole sites as input and -- if solver converges -- provides induced dipoles as output. Since AMOEBA requires the calculations of two sets of induced dipoles generated from two different electric fields (normally called direct (D) and polarization (P)) both electric field and induced dipoles are shaped with an extra dimension and this routine calls the solver twice to solve the two linear systems in the case of AMOEBA FF. Direct electric field and induced dipoles are stored in e(:,:,1)/ipds(:,:,1) while polarization field/dipole are stored in e(:,:,2)/ipds(:,:,2).

Take as input a matrix in which the n-th row stores the atoms that are in the same polarization group as the n-th atom (rows are padded with zeros) and assign to each atom a group index, according to the input information. This is a way of compressing and making more clear the handling of polarization groups.

Perform matrix vector multiplication y = pol*x, where pol is polarizability vector, x and y are column vectors

This subroutine computes the potential generated by the induced point dipoles to a set of arbitrary coordinates, without applying any screening rules. Note: for AMOEBA D dipoles should be used.

This subroutine computes the potential generated by the static multipoles to a set of arbitrary coordinates, without applying any screening rules.

This function allocate and populate array of electrostatic properties of static multipoles at static multipoles sites. It should be called blindly before any calculation that requires V_M2M etc.

Print an array of integers in a well formatted way.

Interface for matrix printing function

TODO Computes the electric potential of a charge $q$ at position $\mathbf{dr}$ from the charge itself. Pre-computed kernel should be provided as input. The result is added to $V$. $$V_q = q \frac{f(r)}{||\mathbf{dr}||}$$

Assign vdw parameters of the QM part from attype and prm file

Computes the missing gradients for QM/MM linkatoms that is bonded terms on QM atoms, LA forces projection on QM and MM atoms. To obtain the correct forces in output, qmg should already contain the QM forces, so that LA forces could be projected on QM and MM force vectors

Allocate a 1-dimensional array of reals

Allocate a 2-dimensional array of reals

Allocate a 3-dimensional array of reals

Free a 1-dimensional array of reals

Free a 2-dimensional array of reals

Free a 3-dimensional array of reals

Reshape a boolean sparse matrix in Yale format to accomodate a larger number of non-zero elements or to trim unused non-zero elements after a guess allocation.

Check which polarizabilities are close enough to 0 to be just excluded from the calculation, and remove them.

Takes as argument an array of group index for each atom, and create a list of atms in each group using the sparse matrix format (saved as Yale format). This is used by cell lists, polarization group etc.

Compute the Jacobian matrix of distance Rij = sqrt((ci(x)-cj(x))2 + (ci(y)-cj(y))2 + (ci(z)-cj(z))**2) Derivatives wrt ci(:) are saved in J_i and wrt cj(:) in J_j; the distance between the two points is also provided in output in Rij.

compute root-mean-square and max norms of a vector.

this routine rotates the atomic multipoles from the molecular frame where they are defined as force field parameters to the lab frame. if required, it also computes the contribution to the forces that stems from the derivatives of the rotation matrices, sometimes referred to as "torques". for the latter task, it uses the field and field gradient from at the multipoles, which is passed in def. for consistency reasons, def is dimensioned (ld_cder,mm_atoms).

given an atom j and the reference atoms jx, jy, and jz, this routine computes the rotation matrix needed to rotate the multipoles on the i-th atom from the molecular frame to the lab frame. if required, it also return the derivative of the rotation matrices with respect to the coordinates of all the atoms involved in its definition.

Utility function used to decide if an interaction between sites i and j should be computed and eventually scaled by a factor. This function is intended to be used in $\mathcalO(n^2)$ code, for linear scaling code lists should be built. This is written to minimize code repetitions, all the screening rules are handled in two possible cases: 1. rules based on adjacency matrix 2. rules based on AMOEBA polarization groups

Set the frozen atoms in the current topology, if the the frozen atoms has already been set, it reinitialize the whole list, without taking into account the content of [[top_obj%frozen]]

This subroutine changes the output file for mmpol to a file defined by filename. A file has already been set, close it before proceed.

Subroutine to initialize the screening parameters

Set the verbosity level for the output, this is a library-level function, that changes the behaviour of several I/O functions. It also enforces min/max verbosity levels (currently no output is 0, while debug output is 3).

Skips n lines while reading an input file

This is a simple -- and unefficient -- routine to sort a vector of integers. It is just used during some output to simplify comparisons with older version of the code.

Inplace equivalent of sort_ivec routine.

Create an identity matrix (boolean sparse, represented in Yale format) of dimension $n$.

Decide if a string starts with a letter or not.

Convert string in input from upper case to lower case and return the lower case string as output.

Remove inline comments from a srting s, comments should begin with the string or character contained in comment_char

Initialize arrays for calculation of stretch-bend cross term potential

Compute the stretch-bend cross term potential.
Those terms are computed according the following formula: $$U_{bond/angle} = (k_i \Delta l_i + k_j \Delta l_j) \Delta \theta_{ij}$$ where $\theta_{ij}$ is the angle delimited by the bond $i$ and $j$.
The force constants $k_i$ and $k_j$ are explicitely defined in the FF, while the equilibrium values are the same as for stretching and bending terms.

Get maximum memory usage since last time push in GB, also make it ready for the next push/pull

This routine compute the product between the diagonal of T matrix with x, and add it to y. The product is simply computed by each element of x for its inverse polarizability.

Perform matrix vector multiplication y = TMat*x, where TMat is polarization matrix (precomputed and stored in memory) and x and y are column vectors

Perform matrix vector multiplication y = [TMat-diag(TMat)]*x, where TMat is polarization matrix (precomputed and stored in memory) and x and y are column vectors

Perform matrix vector multiplication y = TMat*x, where TMat is polarization matrix (precomputed and stored in memory) and x and y are column vectors

This function is used to subsequently break a string into tokens. Tokens separators are any number of spaces.
If just the string is provided, the function returns the position of the first printable character;
If also ib is provided it saves the position of the first printable character after position ib (or ib itself) in ib and return the position of the last printable character before the first space after ib. If ntok is specified instead of a single token, ntok are returned. In case of last token hitten -1 is returned.
To divide a string follow the following scheme:
1. ib = tokenize(s)
2. ie = tokenize(s, ib)
3. tok1 = s(ib:ie)
4a. ib = ib+1
4b. ie = tokenize(s, ib)
5. tok2 = s(ib:ie)

This function is used to subsequently break a string into tokens. Tokens separators are any number of spaces.
If just the string is provided, the function returns the position of the first printable character;
If also ib is provided it saves the position of the first printable character after position ib (or ib itself) in ib and return the position of the last printable character before the first space after ib. If ntok is specified instead of a single token, ntok are returned. In case of last token hitten -1 is returned.
To divide a string follow the following scheme:
1. ib = tokenize(s)
2. ie = tokenize(s, ib)
3. tok1 = s(ib:ie)
4a. ib = ib+1
4b. ie = tokenize(s, ib)
5. tok2 = s(ib:ie)

Computes the Jacobian matrix for torsion angle defined by points $\vec{A}$, $\vec{B}$, $\vec{C}$ and $\vec{D}$ (connected in this order). The angle is defined as follow: [ \vec{U} = (\vec{B} - \vec{C}) \times (\vec{D} - \vec{C}) \ !! \vec{T} = (\vec{B} - \vec{A}) \times (\vec{B} - \vec{C}) \ !! cos(\theta) = \vec{U} \cdot \vec{T}

Compute torsion potential

Initialize torsion potential arrays

Compute torsion potential

Compute torsion potential

Initialize torsion-torsion correction potential arrays

Store in module memory the data describing a new torsion-torsion map

Compute torsion potential

Interface to change the coordinates of the system (eg. during a MD or a geometry optimization). This function clears all the relevant, flags and update the needed quantities. All those operations are needed for a correct functionality of the program therefore coordinates should never be updated without passing from this interface.

Initialize Urey-Bradley potential arrays

Compute the Urey-Bradley potential.
This is basically a virtual bond, with its stretching harminic potential that connect two otherwise un-connected bonds. The same potential formula used for normal stretching is used.

This function is used to know if the library is compiled using integer 8 or 4 bytes long.

Compute the dispersion-repulsion energy using the buffered 7-14 potential. Details can be found in ref: 10.1021/jp027815

Compute the gradient of vdw energy (using the buffered 7-14 potential ref: 10.1021/jp027815) with respect to the distance between the two atoms (Rij).

Compute the dispersion repulsion geometric gradients for the whole system using a double loop algorithm

Compute the dispersion repulsion energy for the whole system using a double loop algorithm

Compute the dispersion repulsion energy for the whole system using a double loop algorithm

Initialize the non-bonded object allocating the parameters vectors

Compute the dispersion repulsion energy for the whole system using a double loop algorithm

Compute the dispersion repulsion energy for the whole system using a double loop algorithm

Compute the dispersion repulsion energy between two systems vdw1 and vdw2 accounting only for the pairs pairs(1,i)--pairs(2,i) and scaling each interaction by s(i).

Remove the VdW interaction from the specified atom the atom will not interact anymore with any other atom

Computes the matrix operator corresponding to a cross product of an input vector $\vec{A}$ of dimension 3. $$[\vec{A}]_\times = \begin{bmatrix} 0 & -\vec{A}_z & \vec{A}_y \ !! \vec{A}_z & 0 & -\vec{A}_x \ !! -\vec{A}_y & \vec{A}_x & 0 \ !! \end{bmatrix}$$

Computes the derivativative matrix of a versor $\hat{A}$ wrt its generator vector $\vec{A}$. $$\hat{A} = \frac{\vec{A}}{||\vec{A}||}$$ [\frac{\partial \hat{A}}{\partial \vec{A}} = \frac{1}{||\vec{A}||^3} (||\vec{A}||^2 \mathbb{1}_3 - A^\dagger A) = \frac{1}{||\vec{A}||^3} \begin{bmatrix} ||\vec{A}||^2 - \vec{A}_x^2 & - \vec{A}_x \vec{A}_y & - \vec{A}_x \vec{A}_z \ !! - \vec{A}_y \vec{A}_x &