mol_md.h File Reference

Header for molecular dynamics. More...

#include "rwmd.h"
#include "rwmolecule.h"
#include "mol_util.h"
#include "utilities.h"

Functions
double	md_leapfrog (Bmolgroup *molgroup, Bmd *md, int max_iter, double velocitylimit)
	Molecular dynamics using the leapfrog integrator. More...
int	md_zero_forces (Bmolgroup *molgroup)
	Zero all atomic forces. More...
double	md_bond_forces (Bmolgroup *molgroup, double Kbond, int wrap)
	Calculates the covalent bond length forces and energy. More...
double	md_angular_forces (Bmolgroup *molgroup, double Kangle, int wrap)
	Calculates the covalent bond angular forces and energy. More...
double	md_nonbonded_forces (Bmolgroup *molgroup, Bmd *md)
	Calculates the non-bonded forces and energy. More...
int	atom_nonbonded_forces (Batom *atom1, Batom *atom2, Bmd *md, Vector3< double > box)
	Calculates the non-bonded forces and energy between two atoms. More...
double	md_point_force (Bmolgroup *molgroup, Vector3< double > point, double Kpoint, double decay)
	Calculates the atomic forces and energy resulting from a single point force. More...

Detailed Description

Header for molecular dynamics.

Author

Bernard Heymann

Date

Created: 20010828

Modified: 20060412

Function Documentation

atom_nonbonded_forces()

int atom_nonbonded_forces	(	Batom *	atom,
		Batom *	atom2,
		Bmd *	md,
		Vector3< double >	box
	)

Calculates the non-bonded forces and energy between two atoms.

Parameters

*atom	central atom.
*atom2	second atom.
*md	molecular dynamics structure.
box	dynamics box.

Returns

double total non-bonded energy.

The energy is defined as a Lennard-Jones term for the Van der Waals 
interactions based on a reference length, ro, and a Coulomb term for 
electrostatic interactions based on the atomic charges, q1 and q2: 
    Enb = Kvdw*((1/12)*(|ro|/|r|)^12 - (1/6)*(|ro|/|r|)^6) + Kelec*q1*q2/|r|
The force is the derivative of the energy:
    Fnb = -Kvdw*((|ro|/|r|)^12 - (|ro|/|r|)^6)*r/|r| + Kelec*q1*q2*r/|r|^3
where r is the distance vector, Kvdw is the Van der Waals energy constant
and Kelec is the electrostatic energy constant.

md_angular_forces()

double md_angular_forces	(	Bmolgroup *	molgroup,
		double	Kangle,
		int	wrap
	)

Calculates the covalent bond angular forces and energy.

Parameters

*molgroup	molecular structure.
Kangle	bond angle energy constant.
wrap	flag to wrap around periodic boundaries.

Returns

double total bond angle energy.

The energy is defined as a harmonic function around the reference 
bond angle, a0:
    Ea = Ka*(cos(a0)-r1*r2/(|r1|*|r2|))^2
The force is the derivative of the energy on the first and last atoms:
    Fa1 = 2*Ka*(cos(a0)-r1*r2/(|r1|*|r2|))/(|r1|*|r2|) * ((r1*r2/|r1|)*r1-r2)
    Fa3 = 2*Ka*(cos(a0)-r1*r2/(|r1|*|r2|))/(|r1|*|r2|) * ((r1*r2/|r2|)*r2-r1)
where r1 is the vector from atom 2 to atom 1, r2 is the vector from
atom 2 to atom 3, and Ka is the bond angle energy constant.

md_bond_forces()

double md_bond_forces	(	Bmolgroup *	molgroup,
		double	Kbond,
		int	wrap
	)

Calculates the covalent bond length forces and energy.

Parameters

*molgroup	molecular structure.
Kbond	bond energy constant.
wrap	flag to wrap around periodic boundaries.

Returns

double total bond length energy.

The energy is defined as a harmonic function around the reference 
bond length, |ro|:
    Eb = Kb*(|r|-|ro|)^2
The force is the derivative of the energy:
    Fb = -2*Kb*(|r|-|ro|)*r/|r|
where r is the distance vector and Kb is the bond energy constant.

md_leapfrog()

double md_leapfrog	(	Bmolgroup *	molgroup,
		Bmd *	md,
		int	max_iter,
		double	velocitylimit
	)

Molecular dynamics using the leapfrog integrator.

Parameters

*molgroup	molecular structure.
*md	molecular dynamics parameters.
max_iter	maximum number of iterations to run.
velocitylimit	limit on velocity per time step.

Returns

double energy.

Leapfrog integration for any coordinate x, velocity vx and force Fx:
    x(t+1) = x(t) + vx(t+1) * dt
    vx(t+1) = (Fx(t) * dt/m + vx(t)) * kf
    where
        kf: friction constant (1=no friction)
        dt: time step
        m: atomic mass
The velocity is limited each time step to damp chaotic oscillations.

md_nonbonded_forces()

double md_nonbonded_forces	(	Bmolgroup *	molgroup,
		Bmd *	md
	)

Calculates the non-bonded forces and energy.

Parameters

*molgroup	molecular structure.
*md	molecular dynamics structure.

Returns

double total non-bonded energy.

The energy is defined as a Lennard-Jones term for the Van der Waals 
interactions based on a reference length, ro, and a Coulomb term for 
electrostatic interactions based on the atomic charges, q1 and q2: 
    Enb = Kvdw*((1/12)*(|ro|/|r|)^12 - (1/6)*(|ro|/|r|)^6) + Kelec*q1*q2/|r|
The force is the derivative of the energy:
    Fnb = -Kvdw*((|ro|/|r|)^12 - (|ro|/|r|)^6)*r/|r| + Kelec*q1*q2*r/|r|^3
where r is the distance vector, Kvdw is the Van der Waals energy constant
and Kelec is the electrostatic energy constant.

md_point_force()

double md_point_force	(	Bmolgroup *	molgroup,
		Vector3< double >	point,
		double	Kpoint,
		double	decay
	)

Calculates the atomic forces and energy resulting from a single point force.

Parameters

*molgroup	molecular structure.
point	center of point force.
Kpoint	point force constant.
decay	energy decay with distance.

Returns

double point force energy.

The energy is defined as an exponential decay over distance from the 
center of the point force:
    Ep = Kp * exp(-decay*dist)
The force is the derivative of the energy:
    Fp = Kp * decay * dir * exp(-decay*dist)
where Kp is the point force constant, dist is the distance of the atom 
from the center of the point force, decay is the energy decay with distance
from the point force center, and dir is the normalized direction vector
pointing from the point force center to the atom, indicating the direction
of force.

md_zero_forces()

int md_zero_forces

(

Bmolgroup *

molgroup

)

Zero all atomic forces.

Parameters

*molgroup

molecular structure.

Returns

int 0.