[brahms-md] r52 committed - Added shifted-force dipole-dipole potential, forces and torques. Impro...
brah...@googlecode.com
Author: orsimario
Date: Wed Jul 20 07:59:24 2011
Log: Added shifted-force dipole-dipole potential, forces and torques.
Improved Lammps interface.
http://code.google.com/p/brahms-md/source/detail?r=52
Modified:
/trunk/src/forceNbdES.c
/trunk/src/forceNbdMain.c
/trunk/src/in_namelist.c
/trunk/src/lammps.c
/trunk/src/main.c
/trunk/src/output.c
/trunk/src/topology.c
=======================================
--- /trunk/src/forceNbdES.c Tue Jun 28 02:10:55 2011
+++ /trunk/src/forceNbdES.c Wed Jul 20 07:59:24 2011
@@ -62,46 +62,52 @@
VScale( *torqueDip, elMagOvrrr * sfFac );
}
-/* Dipolar interaction [allen & tildesley, appendix C] */
-void CalcDipDipInteraction( real *ijEnergy, // out - potential energy
+/* Dipolar interaction [allen & tildesley, appendix C] with shifted-force
truncation */
+void CalcDipDipInteraction_SF( real *ijEnergy, // out - potential energy
VecR *ijForce, // out - force
VecR *ijTorque, // out - torque
VecR *jiTorque, // out - torque
real dipMagProd, // in - product of dipole magnitudes
- const VecR iOrient, // in - "i" orientation (unit) vector
- const VecR jOrient, // in - "j" orientation (unit) vector
+ const VecR ei, // in - "i" orientation (unit) vector
+ const VecR ej, // in - "j" orientation (unit) vector
+ const real rCut, // in - cutoff distance
const VecR rij, // in - intersite distance vector
const real rr ) // in - intersite distance squared
{
- real cosGamma, cosTheta_i, cosTheta_j, r, rrr;
- VecR cosVec1, cosVec2, cosVec, dipTrq1, dipTrq2, rCosVec, rijUnit;
+ real cosGamma, cosTheta_i, cosTheta_j, r, rrr, r_over_rCut,
r_over_rCut_4, fac;
+ VecR cosVec1, cosVec2, cosVec, cosVecFull, e_cross_e, cos_e_cross_r,
rCosVec, rijUnit;
dipMagProd *= ONE_OVER_4_PI_EPSILON0_IN_BRAHMS_UNITS; // scaling
according to the electric conversion factor
r = sqrt ( rr );
rrr = rr * r;
+ r_over_rCut = r/rCut;
+ r_over_rCut_4 = r_over_rCut*r_over_rCut*r_over_rCut*r_over_rCut;
+ fac = 1 - 4*Cube(r_over_rCut) + 3*r_over_rCut_4;
VSCopy( rijUnit, 1./ r, rij );
- cosGamma = VDot( iOrient, jOrient ); /* C.15 */
- cosTheta_i = VDot( iOrient, rijUnit );
- cosTheta_j = VDot( jOrient, rijUnit );
- *ijEnergy = dipMagProd * ( cosGamma - 3. * cosTheta_i * cosTheta_j ) /
rrr; /* (C.16) */
- /* dipolar force calculation - from C.26 a, allen & tildesley */
- VSCopy( cosVec1, cosTheta_j, iOrient );
- VSCopy( cosVec2, cosTheta_i, jOrient );
+ cosGamma = VDot( ei, ej ); /* C.15 */
+ cosTheta_i = VDot( ei, rijUnit );
+ cosTheta_j = VDot( ej, rijUnit );
+ *ijEnergy = dipMagProd * ( cosGamma - 3. * cosTheta_i * cosTheta_j ) *
fac / rrr;
+ /* dipolar force calculation - adapted from C.26 a, allen & tildesley */
+ VSCopy( cosVec1, cosTheta_j, ei );
+ VSCopy( cosVec2, cosTheta_i, ej );
VAdd( cosVec, cosVec1, cosVec2 );
- VSCopy( rCosVec, ( cosGamma - 5. * cosTheta_i * cosTheta_j ), rijUnit );
- VAdd( *ijForce, rCosVec, cosVec );
+ VSAdd( cosVecFull, cosVec, -2*cosTheta_i*cosTheta_j, rijUnit);
+ VScale( cosVecFull, fac );
+ VSCopy( rCosVec, (cosGamma-3*cosTheta_i*cosTheta_j)*(1-r_over_rCut_4),
rijUnit );
+ VAdd( *ijForce, rCosVec, cosVecFull );
VScale( *ijForce, 3. * dipMagProd / Sqr( rr ) );
- /* dipolar torque calculation - from C.26 b-c, allen & tildesley */
- VCross( dipTrq1, iOrient, jOrient );
- VCross( dipTrq2, iOrient, rij );
- VScale( dipTrq2, 3. * cosTheta_j / r );
- VSub( *ijTorque, dipTrq2, dipTrq1 );
- VScale( *ijTorque, dipMagProd / rrr );
- VCross( dipTrq1, jOrient, iOrient );
- VCross( dipTrq2, jOrient, rij ); /* doubt: why not rji? */
- VScale( dipTrq2, 3. * cosTheta_i / r );
- VSub( *jiTorque, dipTrq2, dipTrq1 );
- VScale( *jiTorque, dipMagProd / rrr );
+ /* dipolar torque calculation - adapted from C.26 b-c, allen &
tildesley */
+ VCross( e_cross_e, ei, ej );
+ VCross( cos_e_cross_r, ei, rij );
+ VScale( cos_e_cross_r, 3. * cosTheta_j / r );
+ VSub( *ijTorque, cos_e_cross_r, e_cross_e );
+ VScale( *ijTorque, dipMagProd * fac / rrr );
+ VCross( e_cross_e, ej, ei );
+ VCross( cos_e_cross_r, ej, rij );
+ VScale( cos_e_cross_r, 3. * cosTheta_i / r );
+ VSub( *jiTorque, cos_e_cross_r, e_cross_e );
+ VScale( *jiTorque, dipMagProd * fac / rrr );
}
/* Dipolar interaction [allen & tildesley, appendix C] with cubic
switching for rSwitch < r < rCut */
@@ -110,14 +116,14 @@
VecR *ijTorque, // out - torque
VecR *jiTorque, // out - torque
real dipMagProd, // in - product of dipole magnitudes
- const VecR iOrient, // in - "i" orientation (unit) vector
- const VecR jOrient, // in - "j" orientation (unit) vector
+ const VecR ei, // in - "i" orientation (unit) vector
+ const VecR ej, // in - "j" orientation (unit) vector
const real rCut, // in - cutoff distance
const VecR rij, // in - intersite distance vector
const real rr ) // in - intersite distance squared
{
real cosGamma, cosTheta_i, cosTheta_j, rC_rS_3, r, rrr, rC_r, s, sFactor;
- VecR cosVec1, cosVec2, cosVec, dipTrq1, dipTrq2, rCosVec, rijUnit,
switchGradient;
+ VecR cosVec1, cosVec2, cosVec, e_cross_e, cos_e_cross_r, rCosVec,
rijUnit, switchGradient;
const real rSwitch = 0.9 * rCut; // [Leach, Molecular modelling, p.332,
2nd ed]
@@ -125,27 +131,27 @@
r = sqrt ( rr );
rrr = rr * r;
VSCopy( rijUnit, 1./ r, rij );
- cosGamma = VDot( iOrient, jOrient ); /* C.15 */
- cosTheta_i = VDot( iOrient, rijUnit );
- cosTheta_j = VDot( jOrient, rijUnit );
+ cosGamma = VDot( ei, ej ); /* C.15 */
+ cosTheta_i = VDot( ei, rijUnit );
+ cosTheta_j = VDot( ej, rijUnit );
*ijEnergy = dipMagProd * ( cosGamma - 3. * cosTheta_i * cosTheta_j ) /
rrr; /* (C.16) */
/* dipolar force calculation - from C.26 a, allen & tildesley */
- VSCopy( cosVec1, cosTheta_j, iOrient );
- VSCopy( cosVec2, cosTheta_i, jOrient );
+ VSCopy( cosVec1, cosTheta_j, ei );
+ VSCopy( cosVec2, cosTheta_i, ej );
VAdd( cosVec, cosVec1, cosVec2 );
VSCopy( rCosVec, ( cosGamma - 5. * cosTheta_i * cosTheta_j ), rijUnit );
VAdd( *ijForce, rCosVec, cosVec );
VScale( *ijForce, 3. * dipMagProd / Sqr( rr ) );
/* dipolar torque calculation - from C.26 b-c, allen & tildesley */
- VCross( dipTrq1, iOrient, jOrient );
- VCross( dipTrq2, iOrient, rij );
- VScale( dipTrq2, 3. * cosTheta_j / r );
- VSub( *ijTorque, dipTrq2, dipTrq1 );
+ VCross( e_cross_e, ei, ej );
+ VCross( cos_e_cross_r, ei, rij );
+ VScale( cos_e_cross_r, 3. * cosTheta_j / r );
+ VSub( *ijTorque, cos_e_cross_r, e_cross_e );
VScale( *ijTorque, dipMagProd / rrr );
- VCross( dipTrq1, jOrient, iOrient );
- VCross( dipTrq2, jOrient, rij ); /* doubt: why not rji? */
- VScale( dipTrq2, 3. * cosTheta_i / r );
- VSub( *jiTorque, dipTrq2, dipTrq1 );
+ VCross( e_cross_e, ej, ei );
+ VCross( cos_e_cross_r, ej, rij ); /* doubt: why not rji? */
+ VScale( cos_e_cross_r, 3. * cosTheta_i / r );
+ VSub( *jiTorque, cos_e_cross_r, e_cross_e );
VScale( *jiTorque, dipMagProd / rrr );
if ( r > rSwitch ) { /* switching for rSwitch < r < rCut*/
rC_r = rCut - r;
@@ -160,3 +166,45 @@
*ijEnergy *= s;
}
}
+
+/* Dipolar interaction, straight cutoff [allen & tildesley, appendix C] */
+void CalcDipDipInteraction( real *ijEnergy, // out - potential energy
+ VecR *ijForce, // out - force
+ VecR *ijTorque, // out - torque
+ VecR *jiTorque, // out - torque
+ real dipMagProd, // in - product of dipole magnitudes
+ const VecR ei, // in - "i" orientation (unit) vector
+ const VecR ej, // in - "j" orientation (unit) vector
+ const VecR rij, // in - intersite distance vector
+ const real rr ) // in - intersite distance squared
+{
+ real cosGamma, cosTheta_i, cosTheta_j, r, rrr;
+ VecR cosVec1, cosVec2, cosVec, e_cross_e, cos_e_cross_r, rCosVec,
rijUnit;
+
+ dipMagProd *= ONE_OVER_4_PI_EPSILON0_IN_BRAHMS_UNITS; // scaling
according to the electric conversion factor
+ r = sqrt ( rr );
+ rrr = rr * r;
+ VSCopy( rijUnit, 1./ r, rij );
+ cosGamma = VDot( ei, ej ); /* C.15 */
+ cosTheta_i = VDot( ei, rijUnit );
+ cosTheta_j = VDot( ej, rijUnit );
+ *ijEnergy = dipMagProd * ( cosGamma - 3. * cosTheta_i * cosTheta_j ) /
rrr; /* (C.16) */
+ /* dipolar force calculation - from C.26 a, allen & tildesley */
+ VSCopy( cosVec1, cosTheta_j, ei );
+ VSCopy( cosVec2, cosTheta_i, ej );
+ VAdd( cosVec, cosVec1, cosVec2 );
+ VSCopy( rCosVec, ( cosGamma - 5. * cosTheta_i * cosTheta_j ), rijUnit );
+ VAdd( *ijForce, rCosVec, cosVec );
+ VScale( *ijForce, 3. * dipMagProd / Sqr( rr ) );
+ /* dipolar torque calculation - from C.26 b-c, allen & tildesley */
+ VCross( e_cross_e, ei, ej );
+ VCross( cos_e_cross_r, ei, rij );
+ VScale( cos_e_cross_r, 3. * cosTheta_j / r );
+ VSub( *ijTorque, cos_e_cross_r, e_cross_e );
+ VScale( *ijTorque, dipMagProd / rrr );
+ VCross( e_cross_e, ej, ei );
+ VCross( cos_e_cross_r, ej, rij ); /* doubt: why not rji? */
+ VScale( cos_e_cross_r, 3. * cosTheta_i / r );
+ VSub( *jiTorque, cos_e_cross_r, e_cross_e );
+ VScale( *jiTorque, dipMagProd / rrr );
+}
=======================================
--- /trunk/src/forceNbdMain.c Tue Jun 28 02:10:55 2011
+++ /trunk/src/forceNbdMain.c Wed Jul 20 07:59:24 2011
@@ -28,6 +28,7 @@
void CalcChgChgInteraction_SF( real*, VecR*, real, real, real, VecR );
void CalcDipDipInteraction( real*, VecR*, VecR*, VecR*, real, VecR, VecR,
VecR, real );
void CalcChgDipInteraction_SF( real*, VecR*, VecR*, real, VecR, real,
VecR, real );
+void CalcDipDipInteraction_SF( real*, VecR*, VecR*, VecR*, real, VecR,
VecR, real, VecR, real );
void CalcDipDipInteraction_CS( real*, VecR*, VecR*, VecR*, real, VecR,
VecR, real, VecR, real );
void CalcLennJonesInteraction( real*, VecR*, real, real, real, VecR );
void CalcLennJonesInteraction_SF( real*, VecR*, real, real, real, real,
VecR );
@@ -61,8 +62,8 @@
if ( interactionType[ n ] == WAT_WAT ) { /* water-water interaction
(special case) */
if ( rr < rrCutWatWat ) {
//
CalcDipDipInteraction(&ijPotEn,&ijForceTot,&ijTorqueTot,&jiTorqueTot,watDipoleSq,siteOrientVec[i],siteOrientVec[j],rij,rr);
- CalcDipDipInteraction_CS( &ijPotEn, &ijForceTot, &ijTorqueTot,
&jiTorqueTot, watDipoleSq,
- siteOrientVec[i], siteOrientVec[j], rCutWatWat, rij, rr );
+// CalcDipDipInteraction_CS( &ijPotEn, &ijForceTot, &ijTorqueTot,
&jiTorqueTot, watDipoleSq, siteOrientVec[i], siteOrientVec[j], rCutWatWat,
rij, rr );
+ CalcDipDipInteraction_SF( &ijPotEn, &ijForceTot, &ijTorqueTot,
&jiTorqueTot, watDipoleSq, siteOrientVec[i], siteOrientVec[j], rCutWatWat,
rij, rr );
potEnWatSum += ijPotEn;
// CalcLennJonesInteraction( &ijPotEn, &ijForce, sigWatSq, epsWat,
rr, rij );
CalcLennJonesInteraction_SF( &ijPotEn, &ijForce, sigWatSq, epsWat, rr,
rrCutWatWat, rij );
=======================================
--- /trunk/src/in_namelist.c Thu Jan 6 06:46:39 2011
+++ /trunk/src/in_namelist.c Wed Jul 20 07:59:24 2011
@@ -1,4 +1,4 @@
-/* Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010 Mario Orsi
+/* Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010, 2011 Mario Orsi
This file is part of Brahms.
Brahms is free software: you can redistribute it and/or modify it under
the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or (at
your option) any later version.
@@ -49,8 +49,8 @@
extern VecR region, regionAdjusted, solInertia;
extern VecI initUcell, initHalfCellWat;
extern real scaleRefSpringLength;
-extern real sigChol, sigAmine, sigPhos, sigEst, sigGly, sigWat;
-extern real epsChol, epsAmine, epsPhos, epsGly, epsWat;
+extern real sigChol, sigAmine, sigPhos, sigEst, sigGly, sigWat, sigCharge;
+extern real epsChol, epsAmine, epsPhos, epsGly, epsWat, epsCharge;
extern real hBondAmineGly, hBondWatWat;
extern real hBondWatEst, hBondWatPhos, hBondWatGly;
extern real hBondAmineWat, hBondAminePhos, hBondAmineEst;
@@ -144,12 +144,14 @@
NameR( sigGly ),
NameR( sigEst ),
NameR( sigTail ),
+ NameR( sigCharge ),
NameR( epsChol ),
NameR( epsAmine ),
NameR( epsPhos ),
NameR( epsGly ),
NameR( epsEst ),
NameR( epsTail ),
+ NameR( epsCharge ),
NameR( hBondWatPhos ),
NameR( hBondWatGly ),
NameR( hBondWatEst ),
=======================================
--- /trunk/src/lammps.c Fri Jul 8 09:15:59 2011
+++ /trunk/src/lammps.c Wed Jul 20 07:59:24 2011
@@ -22,7 +22,7 @@
static const int nAtomsPerLipid = 21; // assuming 15-site lipid models
with 2 added atoms for each gly and est site
static const int nBondsPerLipid = 23; // assuming 15-site lipid models: 14
(main) + 6 (tip-tail charges) + 3 (orientation-restraining)
-static const int nBondTypes = 10; // types defined below:
+static const int nBondTypes = 9; // types defined below:
static const int CHOL_PHOS = 1;
static const int PHOS_GLY = 2;
static const int GLY_EST =3;
@@ -339,33 +339,35 @@
}
fprintf(fp, "\n");
- // bond_coeff for 'harmonic' style
- // syntax: bond_coeff bondType rigidityConstant referenceLength
- // NOTE: in LAMMPS, this potential is defined as K(r-r0)^2, so K_LAMMPS
= k_brahms/2
- fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", CHOL_PHOS,
springRigidity_LAMMPS_units/2,
0.9*sig[CHOLINE_TYPE][PHOSPHATE_TYPE]*nm_IN_Angstrom);
- fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", PHOS_GLY,
springRigidity_LAMMPS_units/2,
0.9*sig[PHOSPHATE_TYPE][GLYCEROL_TYPE]*nm_IN_Angstrom);
- fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", GLY_EST,
springRigidity_LAMMPS_units/2,
0.9*sig[GLYCEROL_TYPE][ESTER_TYPE]*nm_IN_Angstrom);
- fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", EST_CHAIN,
springRigidity_LAMMPS_units/2,
0.9*sig[ESTER_TYPE][TAIL_TYPE]*nm_IN_Angstrom);
- fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", CHAIN_CHAIN,
springRigidity_LAMMPS_units/2,
0.9*sig[TAIL_TYPE][TAIL_TYPE]*nm_IN_Angstrom);
- // fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", AMINE_PHOS,
springRigidity_LAMMPS_units/2,
0.9*sig[PHOSPHATE_TYPE][AMINE_TYPE]*nm_IN_Angstrom);
- fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", GLYPLUS_PHOS,
springRigidity_LAMMPS_units/2,
(sig[PHOSPHATE_TYPE][GLYCEROL_TYPE]-sigma[GLYCEROL_TYPE]/4)*nm_IN_Angstrom);
- fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", ESTPLUS_CHAIN,
springRigidity_LAMMPS_units/2,
(sig[ESTER_TYPE][TAIL_TYPE]-sigma[ESTER_TYPE]/4)*nm_IN_Angstrom);
- fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", GLYCHARGE_COM,
springRigidity_LAMMPS_units/2, sigma[GLYCEROL_TYPE]/4*nm_IN_Angstrom);
- fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", ESTCHARGE_COM,
springRigidity_LAMMPS_units/2, sigma[ESTER_TYPE]/4*nm_IN_Angstrom);
- /* fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 10,
springRigidity_LAMMPS_units/2, 0.9*sig[][]*nm_IN_Angstrom); */
- /* fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 11,
springRigidity_LAMMPS_units/2, 0.9*sig[][]*nm_IN_Angstrom); */
- /* fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 12,
springRigidity_LAMMPS_units/2, 0.9*sig[][]*nm_IN_Angstrom); */
- /* fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 13,
springRigidity_LAMMPS_units/2, 0.9*sig[][]*nm_IN_Angstrom); */
- /* fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 14,
springRigidity_LAMMPS_units/2, 0.9*sig[][]*nm_IN_Angstrom); */
- fprintf(fp, "\n");
-
- // angle_coeff for 'cosine/squared' style
- // syntax: angle_coeff angleType rigidityConstant referenceAngle
- // NOTE: in LAMMPS, this potential is defined as
K[cosTheta-cosTheta0]^2, so K_LAMMPS = k_brahms/2
- fprintf(fp, "angle_coeff%5d%8.4f%8.2f\n", _phos_,
angleRigidity_LAMMPS_units/2, refAngle_CholPhosGly*radians_in_degrees);
- fprintf(fp, "angle_coeff%5d%8.4f%8.2f\n", _gly_,
angleRigidity_LAMMPS_units/2, refAngle_PhosGlyEst*radians_in_degrees);
- fprintf(fp, "angle_coeff%5d%8.4f%8.2f\n", _piSat_,
angleRigidity_LAMMPS_units/2, PI*radians_in_degrees);
- fprintf(fp, "angle_coeff%5d%8.4f%8.2f\n", _cisUnsat_,
angleRigidity_LAMMPS_units/2, refAngle_cisUnsat*radians_in_degrees);
- fprintf(fp, "angle_coeff%5d%8.4f%8.2f\n", _piDipCharges_,
angleRigidity_LAMMPS_units/2, PI*radians_in_degrees);
- fprintf(fp, "\n");
-}
+ if (nSites!=nWaters) { // to prevent segmentation fault when simulating
pure water systems
+ // bond_coeff for 'harmonic' style
+ // syntax: bond_coeff bondType rigidityConstant referenceLength
+ // NOTE: in LAMMPS, this potential is defined as K(r-r0)^2, so
K_LAMMPS = k_brahms/2
+ fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", CHOL_PHOS,
springRigidity_LAMMPS_units/2,
0.9*sig[CHOLINE_TYPE][PHOSPHATE_TYPE]*nm_IN_Angstrom);
+ fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", PHOS_GLY,
springRigidity_LAMMPS_units/2,
0.9*sig[PHOSPHATE_TYPE][GLYCEROL_TYPE]*nm_IN_Angstrom);
+ fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", GLY_EST,
springRigidity_LAMMPS_units/2,
0.9*sig[GLYCEROL_TYPE][ESTER_TYPE]*nm_IN_Angstrom);
+ fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", EST_CHAIN,
springRigidity_LAMMPS_units/2,
0.9*sig[ESTER_TYPE][TAIL_TYPE]*nm_IN_Angstrom);
+ fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", CHAIN_CHAIN,
springRigidity_LAMMPS_units/2,
0.9*sig[TAIL_TYPE][TAIL_TYPE]*nm_IN_Angstrom);
+ // fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", AMINE_PHOS,
springRigidity_LAMMPS_units/2,
0.9*sig[PHOSPHATE_TYPE][AMINE_TYPE]*nm_IN_Angstrom);
+ fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", GLYPLUS_PHOS,
springRigidity_LAMMPS_units/2,
(sig[PHOSPHATE_TYPE][GLYCEROL_TYPE]-sigma[GLYCEROL_TYPE]/4)*nm_IN_Angstrom);
+ fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", ESTPLUS_CHAIN,
springRigidity_LAMMPS_units/2,
(sig[ESTER_TYPE][TAIL_TYPE]-sigma[ESTER_TYPE]/4)*nm_IN_Angstrom);
+ fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", GLYCHARGE_COM,
springRigidity_LAMMPS_units/2, sigma[GLYCEROL_TYPE]/4*nm_IN_Angstrom);
+ fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", ESTCHARGE_COM,
springRigidity_LAMMPS_units/2, sigma[ESTER_TYPE]/4*nm_IN_Angstrom);
+ /* fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 10,
springRigidity_LAMMPS_units/2, 0.9*sig[][]*nm_IN_Angstrom); */
+ /* fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 11,
springRigidity_LAMMPS_units/2, 0.9*sig[][]*nm_IN_Angstrom); */
+ /* fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 12,
springRigidity_LAMMPS_units/2, 0.9*sig[][]*nm_IN_Angstrom); */
+ /* fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 13,
springRigidity_LAMMPS_units/2, 0.9*sig[][]*nm_IN_Angstrom); */
+ /* fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 14,
springRigidity_LAMMPS_units/2, 0.9*sig[][]*nm_IN_Angstrom); */
+ fprintf(fp, "\n");
+
+ // angle_coeff for 'cosine/squared' style
+ // syntax: angle_coeff angleType rigidityConstant referenceAngle
+ // NOTE: in LAMMPS, this potential is defined as
K[cosTheta-cosTheta0]^2, so K_LAMMPS = k_brahms/2
+ fprintf(fp, "angle_coeff%5d%8.4f%8.2f\n", _phos_,
angleRigidity_LAMMPS_units/2, refAngle_CholPhosGly*radians_in_degrees);
+ fprintf(fp, "angle_coeff%5d%8.4f%8.2f\n", _gly_,
angleRigidity_LAMMPS_units/2, refAngle_PhosGlyEst*radians_in_degrees);
+ fprintf(fp, "angle_coeff%5d%8.4f%8.2f\n", _piSat_,
angleRigidity_LAMMPS_units/2, PI*radians_in_degrees);
+ fprintf(fp, "angle_coeff%5d%8.4f%8.2f\n", _cisUnsat_,
angleRigidity_LAMMPS_units/2, refAngle_cisUnsat*radians_in_degrees);
+ fprintf(fp, "angle_coeff%5d%8.4f%8.2f\n", _piDipCharges_,
angleRigidity_LAMMPS_units/2, PI*radians_in_degrees);
+ fprintf(fp, "\n");
+ }
+}
=======================================
--- /trunk/src/main.c Wed Jun 29 10:11:43 2011
+++ /trunk/src/main.c Wed Jul 20 07:59:24 2011
@@ -176,7 +176,7 @@
FILE *fPtrAV; // file that records the zConstraintForce, it opened in
setUpJob and closed at the end of main()
/* Topology input parameters - if not supplied in Brahms units, remember
to do the conversion */
-real sigWat, epsWat, sigChol, sigAmine, sigPhos, sigGly, sigEst, epsChol,
epsAmine, epsPhos, epsGly;
+real sigWat, epsWat, sigChol, sigAmine, sigPhos, sigGly, sigEst, epsChol,
epsAmine, epsPhos, epsGly, sigCharge, epsCharge;
real massWat, massChol, massAmine, massPhos, massGly, massEst, massTail;
real inertiaWat, inertiaEst, inertiaGly;
real epsTail, epsEst, estDipole, sigTail;
=======================================
--- /trunk/src/output.c Wed Jul 6 10:04:36 2011
+++ /trunk/src/output.c Wed Jul 20 07:59:24 2011
@@ -34,7 +34,7 @@
void PrintBrahmsVersion()
{
printf("\n**********************\n");
- printf( "* BRAHMS - 06Jul2011 *\n");
+ printf( "* BRAHMS - 20Jul2011 *\n");
printf( "**********************\n\n");
fflush(stdout);
}
=======================================
--- /trunk/src/topology.c Fri Jul 8 09:15:59 2011
+++ /trunk/src/topology.c Wed Jul 20 07:59:24 2011
@@ -1,4 +1,4 @@
-/* Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010 Mario Orsi
+/* Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010, 2011 Mario Orsi
This file is part of Brahms.
Brahms is free software: you can redistribute it and/or modify it under
the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or (at
your option) any later version.
@@ -15,8 +15,8 @@
extern real **sig, **eps, **sigSquared, **productElectroMoments, *mass,
*inertia;
extern real sigWat, epsWat;
extern const real extTemperature, deltaT;
-extern const real sigChol, sigAmine, sigPhos, sigGly, sigEst, sigTail;
-extern const real epsChol, epsAmine, epsPhos, epsEst, epsGly, epsTail;
+extern const real sigChol, sigAmine, sigPhos, sigGly, sigEst, sigTail,
sigCharge;
+extern const real epsChol, epsAmine, epsPhos, epsEst, epsGly, epsTail,
epsCharge;
extern const real massChol, massAmine, massPhos, massGly, soluteMass,
massEst, massTail, massWat;
extern const real inertiaWat, inertiaEst, inertiaGly;
extern const real cholCharge, phosCharge, glyDipole, estDipole, watDipole;
@@ -110,8 +110,8 @@
dipole[ type ] = 0;
break;
case CHARGE_TYPE:
- sigma[ type ] = 0;
- epsilon[ type ] = 0;
+ sigma[ type ] = sigCharge;
+ epsilon[ type ] = epsCharge;
electroMoment[ type ] = 0;
charge[ type ] = 0;
dipole[ type ] = 0;