// This file is part of the ESPResSo distribution (http://www.espresso.mpg.de).
// It is therefore subject to the ESPResSo license agreement which you accepted upon receiving the distribution
// and by which you are legally bound while utilizing this file in any form or way.
// There is NO WARRANTY, not even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// You should have received a copy of that license along with this program;
// if not, refer to http://www.espresso.mpg.de/license.html where its current version can be found, or
// write to Max-Planck-Institute for Polymer Research, Theory Group, PO Box 3148, 55021 Mainz, Germany.
// Copyright (c) 2002-2006; all rights reserved unless otherwise stated.
#ifndef FORCES_H
#define FORCES_H
/** \file forces.h Force calculation.
*
* Responsible:
* Hanjo
*
* \todo Preprocessor switches for all forces (Default: everything is turned on).
* \todo Implement more flexible thermostat, %e.g. which thermostat to use.
*
* For more information see forces.c .
*/
#include
#include "utils.h"
#include "thermostat.h"
#ifdef NPT
#include "pressure.h"
#endif
#include "communication.h"
#ifdef MOLFORCES
#include "topology.h"
#endif
/* include the force files */
#include "p3m.h"
#include "lj.h"
#include "ljgen.h"
#include "steppot.h"
#include "bmhtf-nacl.h"
#include "buckingham.h"
#include "soft_sphere.h"
#include "maggs.h"
#include "tab.h"
#include "ljcos.h"
#include "ljcos2.h"
#include "ljangle.h"
#include "gb.h"
#include "fene.h"
#include "harmonic.h"
#include "subt_lj.h"
#include "subt_gb.h"
#include "angle.h"
#include "dihedral.h"
#include "debye_hueckel.h"
#include "reaction_field.h"
#include "mmm1d.h"
#include "mmm2d.h"
#include "constraint.h"
#include "comforce.h"
#include "comfixed.h"
#include "molforces.h"
#include "morse.h"
#include "elc.h"
#include "adresso.h"
#include "virtual_sites.h"
/** \name Exported Functions */
/************************************************************/
/address@hidden/
/** Calculate forces.
*
* A short list, what the function is doing:
*
* - Initialize forces with: \ref friction_thermo_langevin (ghost forces with zero).
*
- Calculate \ref tcl_bonded "bonded interaction" forces:
* Loop all local particles (not the ghosts).
*
* - FENE
*
- ANGLE (cos bend potential)
*
* - Calculate \ref tcl_non_bonded "non-bonded short range interaction" forces:
* Loop all \ref IA_Neighbor::vList "verlet lists" of all \ref #cells.
*
* - Lennard-Jones.
*
- Buckingham.
*
- Real space part: Coulomb.
*
- Ramp.
*
* - Calculate long range interaction forces:
Uses \ref P3M_calc_kspace_forces.
*
*/
void force_calc();
/** Set forces of all ghosts to zero
*/
void init_forces_ghosts();
/** Calculate non bonded forces between a pair of particles.
@param p1 pointer to particle 1.
@param p2 pointer to particle 2.
@param d vector between p1 and p2.
@param dist distance between p1 and p2.
@param dist2 distance squared between p1 and p2. */
MDINLINE void add_non_bonded_pair_force(Particle *p1, Particle *p2,
double d[3], double dist, double dist2)
{
IA_parameters *ia_params = get_ia_param(p1->p.type,p2->p.type);
double force[3] = { 0., 0., 0. };
double torque1[3] = { 0., 0., 0. };
double torque2[3] = { 0., 0., 0. };
int j;
// printf("Force: %f %f %f \n",p1->f.f[0],p1->f.f[1],p1->f.f[2]);
FORCE_TRACE(fprintf(stderr, "%d: interaction %d<->%d dist %f\n", this_node, p1->p.identity, p2->p.identity, dist));
/***********************************************/
/* thermostat */
/***********************************************/
#ifdef DPD
/* DPD thermostat forces */
if ( thermo_switch & THERMO_DPD ) add_dpd_thermo_pair_force(p1,p2,d,dist,dist2);
#endif
#ifdef INTER_DPD
if ( thermo_switch == THERMO_INTER_DPD ) add_interdpd_pair_force(p1,p2,ia_params,d,dist,dist2);
#endif
/***********************************************/
/* non bonded pair potentials */
/***********************************************/
calc_non_bonded_pair_force(p1,p2,ia_params,d,dist,dist2,force,torque1,torque2);
/***********************************************/
/* short range electrostatics */
/***********************************************/
#ifdef ELECTROSTATICS
if (coulomb.method == COULOMB_DH)
add_dh_coulomb_pair_force(p1,p2,d,dist,force);
if (coulomb.method == COULOMB_RF)
add_rf_coulomb_pair_force(p1,p2,d,dist,force);
#endif
/*********************************************************************/
/* everything before this contributes to the virial pressure in NpT, */
/* but nothing afterwards */
/*********************************************************************/
#ifdef NPT
for (j = 0; j < 3; j++)
if(integ_switch == INTEG_METHOD_NPT_ISO)
nptiso.p_vir[j] += force[j] * d[j];
#endif
/***********************************************/
/* long range electrostatics */
/***********************************************/
#ifdef ELECTROSTATICS
/* real space coulomb */
switch (coulomb.method) {
#ifdef ELP3M
case COULOMB_ELC_P3M: {
add_p3m_coulomb_pair_force(p1->p.q*p2->p.q,d,dist2,dist,force);
// forces from the virtual charges
// they go directly onto the particles, since they are not pairwise forces
if (elc_params.dielectric_contrast_on)
ELC_P3M_dielectric_layers_force_contribution(p1, p2, p1->f.f, p2->f.f);
break;
}
case COULOMB_P3M: {
#ifdef NPT
double eng = add_p3m_coulomb_pair_force(p1->p.q*p2->p.q,d,dist2,dist,force);
if(integ_switch == INTEG_METHOD_NPT_ISO)
nptiso.p_vir[0] += eng;
#else
add_p3m_coulomb_pair_force(p1->p.q*p2->p.q,d,dist2,dist,force);
#endif
break;
}
#endif
case COULOMB_EWALD: {
#ifdef NPT
double eng = add_ewald_coulomb_pair_force(p1,p2,d,dist2,dist,force);
if(integ_switch == INTEG_METHOD_NPT_ISO)
nptiso.p_vir[0] += eng;
#else
add_ewald_coulomb_pair_force(p1,p2,d,dist2,dist,force);
#endif
break;
}
case COULOMB_MMM1D:
add_mmm1d_coulomb_pair_force(p1,p2,d,dist2,dist,force);
break;
case COULOMB_MMM2D:
add_mmm2d_coulomb_pair_force(p1->p.q*p2->p.q,d,dist2,dist,force);
break;
case COULOMB_MAGGS:
if(maggs.yukawa == 1)
add_maggs_yukawa_pair_force(p1,p2,d,dist2,dist,force);
break;
case COULOMB_NONE:
break;
}
#endif /*ifdef ELECTROSTATICS */
/***********************************************/
/* long range magnetostatics */
/***********************************************/
#ifdef MAGNETOSTATICS
/* real space magnetic dipole-dipole */
switch (coulomb.Dmethod) {
#ifdef ELP3M
#ifdef MDLC
case DIPOLAR_MDLC_P3M:
//fall trough
#endif
case DIPOLAR_P3M: {
#ifdef NPT
double eng = add_p3m_dipolar_pair_force(p1,p2,d,dist2,dist,force);
if(integ_switch == INTEG_METHOD_NPT_ISO)
nptiso.p_vir[0] += eng;
#else
add_p3m_dipolar_pair_force(p1,p2,d,dist2,dist,force);
#endif
break;
}
#endif /*ifdef ELP3M */
}
#endif /* ifdef MAGNETOSTATICS */
/***********************************************/
/* add total nonbonded forces to particle */
/***********************************************/
for (j = 0; j < 3; j++) {
p1->f.f[j] += force[j];
p2->f.f[j] -= force[j];
#ifdef ROTATION
p1->f.torque[j] += torque1[j];
p2->f.torque[j] += torque2[j];
#endif
}
// if (p1->p.identity==0 && p2->p.identity == 1) {printf("\n gb01 p1: %f %f %f \n",p1->f.f[0],p1->f.f[1],p1->f.f[2]);}
// if (p1->p.identity==0 && p2->p.identity == 1) {printf(" gb01 p2: %f %f %f \n",p2->f.f[0],p2->f.f[1],p2->f.f[2]);}
// if (p1->p.identity==1 && p2->p.identity == 0) {printf("\n gb10 p1: %f %f %f \n",p1->f.f[0],p1->f.f[1],p1->f.f[2]);}
// if (p1->p.identity==1 && p2->p.identity == 0) {printf(" gb10 p2: %f %f %f \n",p2->f.f[0],p2->f.f[1],p2->f.f[2]);}
}
/** Calculate bonded forces for one particle.
@param p1 particle for which to calculate forces
*/
MDINLINE void add_bonded_force(Particle *p1)
{
// double dx[3], force[3], force2[3], force3[3];
double dx[3]= { 0., 0., 0. };
double force[3] = { 0., 0., 0. };
double force2[3] = { 0., 0., 0. };
double force3[3] = { 0., 0., 0. };
force[0]=force[1]=force[2]=0.0;
#ifdef ROTATION
double torque1[3] = { 0., 0., 0. };
double torque2[3] = { 0., 0., 0. };
#endif
char *errtxt;
Particle *p2, *p3 = NULL, *p4 = NULL;
Bonded_ia_parameters *iaparams;
int i, j, type_num, type, n_partners, bond_broken;
// printf("Force: %f %f %f \n",p1->f.f[0],p1->f.f[1],p1->f.f[2]);
i = 0;
while(ibl.n) {
type_num = p1->bl.e[i++];
iaparams = &bonded_ia_params[type_num];
type = iaparams->type;
n_partners = iaparams->num;
/* fetch particle 2, which is always needed */
p2 = local_particles[p1->bl.e[i++]];
if (!p2) {
// for harmonic spring:
// if cutoff was defined and p2 is not there it is anyway outside the cutoff, see calc_maximal_cutoff()
if ((type==BONDED_IA_HARMONIC)&&(iaparams->p.harmonic.r_cut>0)) return;
errtxt = runtime_error(128 + 2*TCL_INTEGER_SPACE);
ERROR_SPRINTF(errtxt,"{078 bond broken between particles %d and %d (particles not stored on the same node)} ",
p1->p.identity, p1->bl.e[i-1]);
return;
}
/* fetch particle 3 eventually */
if (n_partners >= 2) {
p3 = local_particles[p1->bl.e[i++]];
if (!p3) {
errtxt = runtime_error(128 + 3*TCL_INTEGER_SPACE);
ERROR_SPRINTF(errtxt,"{079 bond broken between particles %d, %d and %d (particles not stored on the same node)} ",
p1->p.identity, p1->bl.e[i-2], p1->bl.e[i-1]);
return;
}
}
/* fetch particle 4 eventually */
if (n_partners >= 3) {
p4 = local_particles[p1->bl.e[i++]];
if (!p4) {
errtxt = runtime_error(128 + 4*TCL_INTEGER_SPACE);
ERROR_SPRINTF(errtxt,"{080 bond broken between particles %d, %d, %d and %d (particles not stored on the same node)} ",
p1->p.identity, p1->bl.e[i-3], p1->bl.e[i-2], p1->bl.e[i-1]);
return;
}
}
if (n_partners == 1) {
/* because of the NPT pressure calculation for pair forces, we need the
1->2 distance vector here. For many body interactions this vector is not needed,
and the pressure calculation not yet clear. */
get_mi_vector(dx, p1->r.p, p2->r.p);
}
switch (type) {
case BONDED_IA_FENE:
bond_broken = calc_fene_pair_force(p1, p2, iaparams, dx, force);
break;
case BONDED_IA_HARMONIC:
bond_broken = calc_harmonic_pair_force(p1, p2, iaparams, dx, force);
break;
#ifdef LENNARD_JONES
case BONDED_IA_SUBT_LJ:
bond_broken = calc_subt_lj_pair_force(p1, p2, iaparams, dx, force);
break;
#endif
#ifdef ROTATION
case BONDED_IA_SUBT_GB:
bond_broken = calc_subt_gb_pair_force(p1, p2, iaparams, dx, force, torque1, torque2);
break;
#endif
#ifdef BOND_ANGLE
case BONDED_IA_ANGLE:
bond_broken = calc_angle_force(p1, p2, p3, iaparams, force, force2);
break;
#endif
case BONDED_IA_DIHEDRAL:
bond_broken = 0;
bond_broken = calc_dihedral_force(p1, p2, p3, p4, iaparams, force, force2, force3);
break;
#ifdef GB_DIHED
case BONDED_IA_GB_DIHEDRAL:
bond_broken = calc_gb_dihedral_torque(p1, p2, iaparams, dx, torque1, torque2,force,force2);
#endif
break;
#ifdef BOND_CONSTRAINT
case BONDED_IA_RIGID_BOND:
//add_rigid_bond_pair_force(p1,p2, iaparams, force, force2);
bond_broken = 0;
force[0]=force[1]=force[2]=0.0;
break;
#endif
#ifdef TABULATED
case BONDED_IA_TABULATED:
switch(iaparams->p.tab.type) {
case TAB_BOND_LENGTH:
bond_broken = calc_tab_bond_force(p1, p2, iaparams, dx, force);
break;
case TAB_BOND_ANGLE:
bond_broken = calc_tab_angle_force(p1, p2, p3, iaparams, force, force2);
break;
case TAB_BOND_DIHEDRAL:
bond_broken = calc_tab_dihedral_force(p1, p2, p3, p4, iaparams, force, force2, force3);
break;
default:
errtxt = runtime_error(128 + TCL_INTEGER_SPACE);
ERROR_SPRINTF(errtxt,"{081 add_bonded_force: tabulated bond type of atom %d unknown\n", p1->p.identity);
return;
}
break;
#endif
#ifdef BOND_VIRTUAL
case BONDED_IA_VIRTUAL_BOND:
bond_broken = 0;
force[0]=force[1]=force[2]=0.0;
break;
#endif
default :
errtxt = runtime_error(128 + TCL_INTEGER_SPACE);
ERROR_SPRINTF(errtxt,"{082 add_bonded_force: bond type of atom %d unknown\n", p1->p.identity);
return;
}
switch (n_partners) {
case 1:
if (bond_broken) {
char *errtext = runtime_error(128 + 2*TCL_INTEGER_SPACE);
ERROR_SPRINTF(errtext,"{083 bond broken between particles %d and %d} ",
p1->p.identity, p2->p.identity);
continue;
}
for (j = 0; j < 3; j++) {
p1->f.f[j] += force[j];
p2->f.f[j] -= force[j];
/* In subt_gb, need to subtract torques */
#ifdef ROTATION
p1->f.torque[j] += torque1[j];
p2->f.torque[j] += torque2[j];
#endif
#ifdef NPT
if(integ_switch == INTEG_METHOD_NPT_ISO)
nptiso.p_vir[j] += force[j] * dx[j];
#endif
}
break;
case 2:
if (bond_broken) {
char *errtext = runtime_error(128 + 3*TCL_INTEGER_SPACE);
ERROR_SPRINTF(errtext,"{084 bond broken between particles %d, %d and %d} ",
p1->p.identity, p2->p.identity, p3->p.identity);
continue;
}
for (j = 0; j < 3; j++) {
p1->f.f[j] += force[j];
p2->f.f[j] += force2[j];
p3->f.f[j] -= (force[j] + force2[j]);
}
break;
case 3:
if (bond_broken) {
char *errtext = runtime_error(128 + 4*TCL_INTEGER_SPACE);
ERROR_SPRINTF(errtext,"{085 bond broken between particles %d, %d, %d and %d} ",
p1->p.identity, p2->p.identity, p3->p.identity, p4->p.identity);
continue;
}
for (j = 0; j < 3; j++) {
p1->f.f[j] += force[j];
p2->f.f[j] += force2[j];
p3->f.f[j] += force3[j];
p4->f.f[j] -= (force[j] + force2[j] + force3[j]);
}
break;
}
}
}
/** add force to another. This is used when collecting ghost forces. */
MDINLINE void add_force(ParticleForce *F_to, ParticleForce *F_add)
{
int i;
for (i = 0; i < 3; i++)
F_to->f[i] += F_add->f[i];
#ifdef ROTATION
for (i = 0; i < 3; i++)
F_to->torque[i] += F_add->torque[i];
#endif
}
/address@hidden/
#endif