[brahms-md] r46 committed - Added 'CalcDipDipInteraction' (dipolar interactions with straight cuto...
brah...@googlecode.com
Author: orsimario
Date: Tue Jun 28 02:10:55 2011
Log: Added 'CalcDipDipInteraction' (dipolar interactions with straight
cutoff) for debugging purposes and comparison to other codes. Renamed a few
variables in 'forceNbdES.c' for improved clarity. Moved obsolete
module 'forceNbdStky.c' into depot.
http://code.google.com/p/brahms-md/source/detail?r=46
Added:
/trunk/src/depot/forceNbdStky.c
Deleted:
/trunk/src/forceNbdStky.c
Modified:
/trunk/src/forceNbdES.c
/trunk/src/forceNbdMain.c
/trunk/src/output.c
/trunk/src/startUp.c
/trunk/src/sysGen.c
=======================================
--- /dev/null
+++ /trunk/src/depot/forceNbdStky.c Tue Jun 28 02:10:55 2011
@@ -0,0 +1,180 @@
+/* Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010 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.
+ Brahms is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.
+ You should have received a copy of the GNU General Public License along
with Brahms. If not, see <http://www.gnu.org/licenses/>. */
+
+// forceNbdStky.c -- functions for evaluation of the SSD "sticky" potential
+
+#include "definitions.h"
+
+/* local (to this module) variable -- values from Fennel and Gezelter, J
Chem Phys, 120, 9175 (2004) */
+static const real v0 = 3.90 * cal_IN_J, // [kJ/mol]
+ w0 = 0.07715;
+static const real halfv0 = 0.5 * 3.90 * cal_IN_J; // v0 / 2 [kJ/mol]
+static const real rL = 0.24, // [nm]
+ rL1 = 0.275; // [nm]
+const real rU = 0.38, // [nm]
+ rU1 = 0.335; // [nm]
+
+/*************************************************************************************************************************************
+ * CalcStickyInteraction -- calculates SSD-SSD "sticky" pair interaction.
Reference: Chandra & Ichiye, J Chem Phys, 11, 2701 (1999). *
+
*************************************************************************************************************************************/
+void CalcStickyInteraction( real *uij, // out - potential energy
+ VecR *stickyForce, // out - force
+ VecR *stickyTorqueij, // out - torque
+ VecR *stickyTorqueji, // out - torque
+ const VecR dr, // in - distance
+ const RMat rMati,
+ const RMat rMatj,
+ const real rr ) // in - r*r
+{
+ real s, s1; /* switching factors */
+ real r, sFactor, rUmr, CrUmrL;
+ real rrr; /* cube( pair distance ) */
+ real wij, wji, w; /* tetrahedral sticky potential energy */
+ real w1ij, w1ji, w1 = 0.; /* tetrahedral sticky potential correction
energy */
+ real fwijFactor, fwjiFactor, rbijzSq, rbjizSq, xxijMinusyyij, /* xij^2 -
yij^2 */
+ xxijMinusyyijDivrrr, /* (xij^2 - yij^2)/rrr */
+ zijDivrr, /* z/rr */
+ zijDivrrr, /* z/rrr */
+ zr6ij, /* zij/rbij - 0.6 */
+ zr8ij, /* zij/rbij + 0.8 */
+ zr62ij, /* (zij/rbij - 0.6)^2 */
+ zr82ij; /* (zij/rbij + 0.8)^2 */
+ real xxjiMinusyyji, /* xji^2 - yji^2 */
+ xxjiMinusyyjiDivrrr, /* (xji^2 - yji^2)/rrr */
+ zjiDivrr, /* z/rr */
+ zjiDivrrr, /* z/rrr */
+ zr6ji, /* zji/rbji - 0.6 */
+ zr8ji, /* zji/rbji + 0.8 */
+ zr62ji, /* (zji/rbji - 0.6)^2 */
+ zr82ji; /* (zji/rbji + 0.8)^2 */
+ real ri, stickyScalingFactorij, stickyScalingFactorji;
+ VecR drji, rbij, rbji, fibwij, fibw1ij, fibwji, fibw1ji, Fib_ij, Fib_ji,
+ FiSij, FiSji, Tibwij, Tibwji, Tibw1ij, Tibw1ji, Ti_ij, Tj_ji,
switchGradient;
+
+ VZero( fibw1ij );
+ VZero( Tibw1ij );
+ VZero( fibw1ji );
+ VZero( Tibw1ji );
+ // rbij is defined as the position of molecule j in the frame fixed on
molecule i, so rbij = rMati * (-drij)
+ VSCopy( drji, -1., dr );
+ MVMul( rbij, rMati.u, drji ); /* rij = Ri * drji */
+ // printf("*** rbij = % f % f % f\n", rbij.x, rbij.y, rbij.z);
+ xxijMinusyyij = Sqr( rbij.x ) - Sqr( rbij.y );
+ r = sqrt ( rr );
+ rrr = r * rr;
+ xxijMinusyyijDivrrr = xxijMinusyyij / rrr;
+ zr6ij = rbij.z / r - 0.6;
+ zr8ij = rbij.z / r + 0.8;
+ zr62ij = Sqr( zr6ij );
+ zr82ij = Sqr( zr8ij );
+ zijDivrr = rbij.z / rr;
+ zijDivrrr = zijDivrr / r;
+ rbijzSq = Sqr( rbij.z );
+ fwijFactor = - 6. * zijDivrr * xxijMinusyyijDivrrr;
+ wij = 2. * xxijMinusyyijDivrrr * rbij.z; /* eq. (8) from Chandra &
Ichiye, J Chem Phys, 11, 2701 (1999) */
+ /* evaluate wij-part of sticky force on i */
+ fibwij.x = 4. * rbij.x * zijDivrrr + fwijFactor * rbij.x; /* (12)
from Chandra & Ichiye, J Chem Phys, 11, 2701 (1999) */
+ fibwij.y = - 4. * rbij.y * zijDivrrr + fwijFactor * rbij.y; /*
(13)[ditto for this and following numbers] */
+ fibwij.z = 2. * xxijMinusyyijDivrrr + fwijFactor * rbij.z; /* (14) */
+ /* evaluate wij-part of sticky torque on i */
+ Tibwij.x = 4. * ( rbij.y * rbijzSq + rbij.y * xxijMinusyyij / 2. ); /*
(21) */
+ Tibwij.y = 4. * ( rbij.x * rbijzSq - rbij.x * xxijMinusyyij / 2. ); /*
(22) */
+ Tibwij.z = - 8. * VProd( rbij ); /*
(23) */
+ VScale( Tibwij, 1. / rrr );
+ MVMul( rbji, rMatj.u, dr ); /* rbji contains dr_ji's coordinates in
molecule-j's frame */
+ xxjiMinusyyji = Sqr( rbji.x ) - Sqr( rbji.y );
+ xxjiMinusyyjiDivrrr = xxjiMinusyyji / rrr;
+ zr6ji = rbji.z / r - 0.6;
+ zr8ji = rbji.z / r + 0.8;
+ zr62ji = Sqr( zr6ji );
+ zr82ji = Sqr( zr8ji );
+ zjiDivrr = rbji.z / rr;
+ zjiDivrrr = zjiDivrr / r;
+ rbjizSq = Sqr( rbji.z );
+ fwjiFactor = - 6. * zjiDivrr * xxjiMinusyyjiDivrrr;
+ wji = 2. * xxjiMinusyyjiDivrrr * rbji.z;
+ w = wij + wji ;
+ /* evaluate wji-part of sticky force on i */
+ fibwji.x = 4. * rbji.x * zjiDivrrr + fwjiFactor * rbji.x; /* (12) */
+ fibwji.y = - 4. * rbji.y * zjiDivrrr + fwjiFactor * rbji.y; /* (13) */
+ fibwji.z = 2. * xxjiMinusyyjiDivrrr + fwjiFactor * rbji.z; /* (14) */
+ /* evaluate wji-part of sitcky torque on j */
+ Tibwji.x = 4. * ( rbji.y * rbjizSq + rbji.y * xxjiMinusyyji / 2.
); /* (21) */
+ Tibwji.y = 4. * ( rbji.x * rbjizSq - rbji.x * xxjiMinusyyji / 2.
); /* (22) */
+ Tibwji.z = - 8. * VProd( rbji
); /* (23) */
+ VScale( Tibwji, 1. / rrr );
+ if ( r > rL ) { /* switching */
+ rUmr = rU - r;
+ CrUmrL = Cube( rU - rL );
+ s = Sqr( rUmr ) * ( rU + 2. * r - 3. * rL ) / CrUmrL;
+ w *= s;
+ VScale( fibwij, s );
+ VScale( fibwji, s );
+ VScale( Tibwij, s );
+ VScale( Tibwji, s );
+ sFactor = 6. * rUmr * ( rL - r ) / ( CrUmrL * r );
+ VSCopy( switchGradient, sFactor, rbij );
+ VVSAdd( fibwij, wij, switchGradient );
+ VSCopy( switchGradient, sFactor, rbji );
+ VVSAdd( fibwji, wji, switchGradient );
+ }
+ if ( r < rU1 ) { /* compute "w1"-part */
+ /* w1ij */
+ ri = 1. / r;
+ w1ij = zr62ij * zr82ij - w0;
+ fibw1ij.x = -2. * rbij.x * zijDivrrr; /* (15) */
+ fibw1ij.y = -2. * rbij.y * zijDivrrr; /* (16) */
+ fibw1ij.z = 2. * ( ri - zijDivrrr * rbij.z ); /* (17) */
+ stickyScalingFactorij = zr62ij * zr8ij + zr6ij * zr82ij;
+ VScale( fibw1ij, stickyScalingFactorij );
+ Tibw1ij.x = 2. * rbij.y * ri; /* (24) */
+ Tibw1ij.y = - 2. * rbij.x * ri; /* (25) */
+ Tibw1ij.z = 0.; /* (26) */
+ VScale( Tibw1ij, stickyScalingFactorij );
+ /* w1ji */
+ w1ji = zr62ji * zr82ji - w0;
+ w1 = w1ij + w1ji;
+ fibw1ji.x = -2. * rbji.x * zjiDivrrr; /* (15) */
+ fibw1ji.y = -2. * rbji.y * zjiDivrrr; /* (16) */
+ fibw1ji.z = 2. * ( ri - zjiDivrrr * rbji.z ); /* (17) */
+ stickyScalingFactorji = zr62ji * zr8ji + zr6ji * zr82ji;
+ VScale( fibw1ji, stickyScalingFactorji );
+ Tibw1ji.x = 2. * rbji.y * ri; /* (24) */
+ Tibw1ji.y = - 2. * rbji.x * ri; /* (25) */
+ Tibw1ji.z = 0.; /* (26) */
+ VScale( Tibw1ji, stickyScalingFactorji );
+ if ( r > rL1 ) { /* switching, for rL1 < r < rU1 */
+ rUmr = rU1 - r;
+ CrUmrL = Cube( rU1 - rL1 );
+ s1 = Sqr( rUmr ) * ( rU1 + 2. * r - 3. * rL1 ) / CrUmrL;
+ w1 *= s1;
+ VScale( fibw1ij, s1 );
+ VScale( fibw1ji, s1 );
+ VScale( Tibw1ij, s1 );
+ VScale( Tibw1ji, s1 );
+ sFactor = 6. * rUmr * ( rL1 - r ) / ( CrUmrL * r );
+ VSCopy( switchGradient, sFactor, rbij );
+ VVSAdd( fibw1ij, w1ij, switchGradient );
+ VSCopy( switchGradient, sFactor, rbji );
+ VVSAdd( fibw1ji, w1ji, switchGradient );
+ }
+ } /* end if ( r < ru1 ) */
+ *uij = halfv0 * ( w + w1 );
+ VAdd( Fib_ij, fibwij, fibw1ij ); /* i-frame */
+ VAdd( Ti_ij, Tibwij, Tibw1ij ); /* i-frame */
+ VAdd( Fib_ji, fibwji, fibw1ji ); /* j-frame */
+ VAdd( Tj_ji, Tibwji, Tibw1ji ); /* j-frame */
+ MVMulT( FiSij, rMati.u, Fib_ij ); /* FiSij = Ri^T * Fib_ij */
+ MVMulT( FiSji, rMatj.u, Fib_ji ); /* FiSji = Rj^T * Fib_ji */
+ VSub( *stickyForce, FiSij, FiSji );
+ MVMulT( *stickyTorqueij, rMati.u, Ti_ij ); /* *stickyTorqueij = Ri^T *
Ti_ij */
+ MVMulT( *stickyTorqueji, rMatj.u, Tj_ji ); /* *stickyTorqueji = Rj^T *
Tj_ji */
+ VScale( *stickyForce, halfv0 ); // 21/11/2006: fixed sign mistake
+ VScale( *stickyTorqueij, halfv0 );
+ VScale( *stickyTorqueji, halfv0 );
+}
=======================================
--- /trunk/src/forceNbdStky.c Thu Jan 6 06:46:39 2011
+++ /dev/null
@@ -1,180 +0,0 @@
-/* Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010 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.
- Brahms is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.
- You should have received a copy of the GNU General Public License along
with Brahms. If not, see <http://www.gnu.org/licenses/>. */
-
-// forceNbdStky.c -- functions for evaluation of the SSD "sticky" potential
-
-#include "definitions.h"
-
-/* local (to this module) variable -- values from Fennel and Gezelter, J
Chem Phys, 120, 9175 (2004) */
-static const real v0 = 3.90 * cal_IN_J, // [kJ/mol]
- w0 = 0.07715;
-static const real halfv0 = 0.5 * 3.90 * cal_IN_J; // v0 / 2 [kJ/mol]
-static const real rL = 0.24, // [nm]
- rL1 = 0.275; // [nm]
-const real rU = 0.38, // [nm]
- rU1 = 0.335; // [nm]
-
-/*************************************************************************************************************************************
- * CalcStickyInteraction -- calculates SSD-SSD "sticky" pair interaction.
Reference: Chandra & Ichiye, J Chem Phys, 11, 2701 (1999). *
-
*************************************************************************************************************************************/
-void CalcStickyInteraction( real *uij, // out - potential energy
- VecR *stickyForce, // out - force
- VecR *stickyTorqueij, // out - torque
- VecR *stickyTorqueji, // out - torque
- const VecR dr, // in - distance
- const RMat rMati,
- const RMat rMatj,
- const real rr ) // in - r*r
-{
- real s, s1; /* switching factors */
- real r, sFactor, rUmr, CrUmrL;
- real rrr; /* cube( pair distance ) */
- real wij, wji, w; /* tetrahedral sticky potential energy */
- real w1ij, w1ji, w1 = 0.; /* tetrahedral sticky potential correction
energy */
- real fwijFactor, fwjiFactor, rbijzSq, rbjizSq, xxijMinusyyij, /* xij^2 -
yij^2 */
- xxijMinusyyijDivrrr, /* (xij^2 - yij^2)/rrr */
- zijDivrr, /* z/rr */
- zijDivrrr, /* z/rrr */
- zr6ij, /* zij/rbij - 0.6 */
- zr8ij, /* zij/rbij + 0.8 */
- zr62ij, /* (zij/rbij - 0.6)^2 */
- zr82ij; /* (zij/rbij + 0.8)^2 */
- real xxjiMinusyyji, /* xji^2 - yji^2 */
- xxjiMinusyyjiDivrrr, /* (xji^2 - yji^2)/rrr */
- zjiDivrr, /* z/rr */
- zjiDivrrr, /* z/rrr */
- zr6ji, /* zji/rbji - 0.6 */
- zr8ji, /* zji/rbji + 0.8 */
- zr62ji, /* (zji/rbji - 0.6)^2 */
- zr82ji; /* (zji/rbji + 0.8)^2 */
- real ri, stickyScalingFactorij, stickyScalingFactorji;
- VecR drji, rbij, rbji, fibwij, fibw1ij, fibwji, fibw1ji, Fib_ij, Fib_ji,
- FiSij, FiSji, Tibwij, Tibwji, Tibw1ij, Tibw1ji, Ti_ij, Tj_ji,
switchGradient;
-
- VZero( fibw1ij );
- VZero( Tibw1ij );
- VZero( fibw1ji );
- VZero( Tibw1ji );
- // rbij is defined as the position of molecule j in the frame fixed on
molecule i, so rbij = rMati * (-drij)
- VSCopy( drji, -1., dr );
- MVMul( rbij, rMati.u, drji ); /* rij = Ri * drji */
- // printf("*** rbij = % f % f % f\n", rbij.x, rbij.y, rbij.z);
- xxijMinusyyij = Sqr( rbij.x ) - Sqr( rbij.y );
- r = sqrt ( rr );
- rrr = r * rr;
- xxijMinusyyijDivrrr = xxijMinusyyij / rrr;
- zr6ij = rbij.z / r - 0.6;
- zr8ij = rbij.z / r + 0.8;
- zr62ij = Sqr( zr6ij );
- zr82ij = Sqr( zr8ij );
- zijDivrr = rbij.z / rr;
- zijDivrrr = zijDivrr / r;
- rbijzSq = Sqr( rbij.z );
- fwijFactor = - 6. * zijDivrr * xxijMinusyyijDivrrr;
- wij = 2. * xxijMinusyyijDivrrr * rbij.z; /* eq. (8) from Chandra &
Ichiye, J Chem Phys, 11, 2701 (1999) */
- /* evaluate wij-part of sticky force on i */
- fibwij.x = 4. * rbij.x * zijDivrrr + fwijFactor * rbij.x; /* (12)
from Chandra & Ichiye, J Chem Phys, 11, 2701 (1999) */
- fibwij.y = - 4. * rbij.y * zijDivrrr + fwijFactor * rbij.y; /*
(13)[ditto for this and following numbers] */
- fibwij.z = 2. * xxijMinusyyijDivrrr + fwijFactor * rbij.z; /* (14) */
- /* evaluate wij-part of sticky torque on i */
- Tibwij.x = 4. * ( rbij.y * rbijzSq + rbij.y * xxijMinusyyij / 2. ); /*
(21) */
- Tibwij.y = 4. * ( rbij.x * rbijzSq - rbij.x * xxijMinusyyij / 2. ); /*
(22) */
- Tibwij.z = - 8. * VProd( rbij ); /*
(23) */
- VScale( Tibwij, 1. / rrr );
- MVMul( rbji, rMatj.u, dr ); /* rbji contains dr_ji's coordinates in
molecule-j's frame */
- xxjiMinusyyji = Sqr( rbji.x ) - Sqr( rbji.y );
- xxjiMinusyyjiDivrrr = xxjiMinusyyji / rrr;
- zr6ji = rbji.z / r - 0.6;
- zr8ji = rbji.z / r + 0.8;
- zr62ji = Sqr( zr6ji );
- zr82ji = Sqr( zr8ji );
- zjiDivrr = rbji.z / rr;
- zjiDivrrr = zjiDivrr / r;
- rbjizSq = Sqr( rbji.z );
- fwjiFactor = - 6. * zjiDivrr * xxjiMinusyyjiDivrrr;
- wji = 2. * xxjiMinusyyjiDivrrr * rbji.z;
- w = wij + wji ;
- /* evaluate wji-part of sticky force on i */
- fibwji.x = 4. * rbji.x * zjiDivrrr + fwjiFactor * rbji.x; /* (12) */
- fibwji.y = - 4. * rbji.y * zjiDivrrr + fwjiFactor * rbji.y; /* (13) */
- fibwji.z = 2. * xxjiMinusyyjiDivrrr + fwjiFactor * rbji.z; /* (14) */
- /* evaluate wji-part of sitcky torque on j */
- Tibwji.x = 4. * ( rbji.y * rbjizSq + rbji.y * xxjiMinusyyji / 2.
); /* (21) */
- Tibwji.y = 4. * ( rbji.x * rbjizSq - rbji.x * xxjiMinusyyji / 2.
); /* (22) */
- Tibwji.z = - 8. * VProd( rbji
); /* (23) */
- VScale( Tibwji, 1. / rrr );
- if ( r > rL ) { /* switching */
- rUmr = rU - r;
- CrUmrL = Cube( rU - rL );
- s = Sqr( rUmr ) * ( rU + 2. * r - 3. * rL ) / CrUmrL;
- w *= s;
- VScale( fibwij, s );
- VScale( fibwji, s );
- VScale( Tibwij, s );
- VScale( Tibwji, s );
- sFactor = 6. * rUmr * ( rL - r ) / ( CrUmrL * r );
- VSCopy( switchGradient, sFactor, rbij );
- VVSAdd( fibwij, wij, switchGradient );
- VSCopy( switchGradient, sFactor, rbji );
- VVSAdd( fibwji, wji, switchGradient );
- }
- if ( r < rU1 ) { /* compute "w1"-part */
- /* w1ij */
- ri = 1. / r;
- w1ij = zr62ij * zr82ij - w0;
- fibw1ij.x = -2. * rbij.x * zijDivrrr; /* (15) */
- fibw1ij.y = -2. * rbij.y * zijDivrrr; /* (16) */
- fibw1ij.z = 2. * ( ri - zijDivrrr * rbij.z ); /* (17) */
- stickyScalingFactorij = zr62ij * zr8ij + zr6ij * zr82ij;
- VScale( fibw1ij, stickyScalingFactorij );
- Tibw1ij.x = 2. * rbij.y * ri; /* (24) */
- Tibw1ij.y = - 2. * rbij.x * ri; /* (25) */
- Tibw1ij.z = 0.; /* (26) */
- VScale( Tibw1ij, stickyScalingFactorij );
- /* w1ji */
- w1ji = zr62ji * zr82ji - w0;
- w1 = w1ij + w1ji;
- fibw1ji.x = -2. * rbji.x * zjiDivrrr; /* (15) */
- fibw1ji.y = -2. * rbji.y * zjiDivrrr; /* (16) */
- fibw1ji.z = 2. * ( ri - zjiDivrrr * rbji.z ); /* (17) */
- stickyScalingFactorji = zr62ji * zr8ji + zr6ji * zr82ji;
- VScale( fibw1ji, stickyScalingFactorji );
- Tibw1ji.x = 2. * rbji.y * ri; /* (24) */
- Tibw1ji.y = - 2. * rbji.x * ri; /* (25) */
- Tibw1ji.z = 0.; /* (26) */
- VScale( Tibw1ji, stickyScalingFactorji );
- if ( r > rL1 ) { /* switching, for rL1 < r < rU1 */
- rUmr = rU1 - r;
- CrUmrL = Cube( rU1 - rL1 );
- s1 = Sqr( rUmr ) * ( rU1 + 2. * r - 3. * rL1 ) / CrUmrL;
- w1 *= s1;
- VScale( fibw1ij, s1 );
- VScale( fibw1ji, s1 );
- VScale( Tibw1ij, s1 );
- VScale( Tibw1ji, s1 );
- sFactor = 6. * rUmr * ( rL1 - r ) / ( CrUmrL * r );
- VSCopy( switchGradient, sFactor, rbij );
- VVSAdd( fibw1ij, w1ij, switchGradient );
- VSCopy( switchGradient, sFactor, rbji );
- VVSAdd( fibw1ji, w1ji, switchGradient );
- }
- } /* end if ( r < ru1 ) */
- *uij = halfv0 * ( w + w1 );
- VAdd( Fib_ij, fibwij, fibw1ij ); /* i-frame */
- VAdd( Ti_ij, Tibwij, Tibw1ij ); /* i-frame */
- VAdd( Fib_ji, fibwji, fibw1ji ); /* j-frame */
- VAdd( Tj_ji, Tibwji, Tibw1ji ); /* j-frame */
- MVMulT( FiSij, rMati.u, Fib_ij ); /* FiSij = Ri^T * Fib_ij */
- MVMulT( FiSji, rMatj.u, Fib_ji ); /* FiSji = Rj^T * Fib_ji */
- VSub( *stickyForce, FiSij, FiSji );
- MVMulT( *stickyTorqueij, rMati.u, Ti_ij ); /* *stickyTorqueij = Ri^T *
Ti_ij */
- MVMulT( *stickyTorqueji, rMatj.u, Tj_ji ); /* *stickyTorqueji = Rj^T *
Tj_ji */
- VScale( *stickyForce, halfv0 ); // 21/11/2006: fixed sign mistake
- VScale( *stickyTorqueij, halfv0 );
- VScale( *stickyTorqueji, halfv0 );
-}
=======================================
--- /trunk/src/forceNbdES.c Thu Jan 6 06:46:39 2011
+++ /trunk/src/forceNbdES.c Tue Jun 28 02:10:55 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.
@@ -61,6 +61,48 @@
VCross( *torqueDip, rChgDip, e );
VScale( *torqueDip, elMagOvrrr * sfFac );
}
+
+/* Dipolar interaction [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 iOrient, // in - "i" orientation (unit) vector
+ const VecR jOrient, // 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, dipTrq1, dipTrq2, 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( 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 );
+ 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 );
+ 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 interaction [allen & tildesley, appendix C] with cubic
switching for rSwitch < r < rCut */
void CalcDipDipInteraction_CS( real *ijEnergy, // out - potential energy
@@ -74,7 +116,7 @@
const VecR rij, // in - intersite distance vector
const real rr ) // in - intersite distance squared
{
- real cosGamma, cosTheta_i, cosTheta_j, CrUmrL, r, rrr, rUmr, s, sFactor;
+ real cosGamma, cosTheta_i, cosTheta_j, rC_rS_3, r, rrr, rC_r, s, sFactor;
VecR cosVec1, cosVec2, cosVec, dipTrq1, dipTrq2, rCosVec, rijUnit,
switchGradient;
const real rSwitch = 0.9 * rCut; // [Leach, Molecular modelling, p.332,
2nd ed]
@@ -106,11 +148,11 @@
VSub( *jiTorque, dipTrq2, dipTrq1 );
VScale( *jiTorque, dipMagProd / rrr );
if ( r > rSwitch ) { /* switching for rSwitch < r < rCut*/
- rUmr = rCut - r;
- CrUmrL = Cube( rCut - rSwitch );
- sFactor = 6. * rUmr * ( rSwitch - r ) / ( CrUmrL * r );
+ rC_r = rCut - r;
+ rC_rS_3 = Cube( rCut - rSwitch );
+ sFactor = 6. * rC_r * ( rSwitch - r ) / ( rC_rS_3 * r );
VSCopy( switchGradient, sFactor, rij );
- s = Sqr( rUmr ) * ( rCut + 2. * r - 3. * rSwitch ) / CrUmrL;
+ s = Sqr( rC_r ) * ( rCut + 2. * r - 3. * rSwitch ) / rC_rS_3;
VScale( *ijForce, s );
VScale( *ijTorque, s );
VScale( *jiTorque, s );
=======================================
--- /trunk/src/forceNbdMain.c Thu Jan 6 06:46:39 2011
+++ /trunk/src/forceNbdMain.c Tue Jun 28 02:10:55 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.
@@ -26,6 +26,7 @@
/* external function declarations: */
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_CS( real*, VecR*, VecR*, VecR*, real, VecR,
VecR, real, VecR, real );
void CalcLennJonesInteraction( real*, VecR*, real, real, real, VecR );
@@ -59,6 +60,7 @@
if ( rr < rrCutMax ) {
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 );
potEnWatSum += ijPotEn;
=======================================
--- /trunk/src/output.c Wed Feb 9 07:48:13 2011
+++ /trunk/src/output.c Tue Jun 28 02:10:55 2011
@@ -34,7 +34,7 @@
void PrintBrahmsVersion()
{
printf("\n**********************\n");
- printf( "* BRAHMS - 09Feb2011 *\n");
+ printf( "* BRAHMS - 28Jun2011 *\n");
printf( "**********************\n\n");
fflush(stdout);
}
=======================================
--- /trunk/src/startUp.c Wed Feb 9 07:48:13 2011
+++ /trunk/src/startUp.c Tue Jun 28 02:10:55 2011
@@ -18,7 +18,7 @@
extern Prop boxVolume, xyArea;
extern Site *site; /* array of Site, represents the state of the system */
extern Atom *atom; /* array of Atom, represents solute atoms */
-extern real *mass, *inertia, *rLipid, timeNow, soluteMass, totalMass;;
+extern real *mass, *inertia, *rLipid, timeNow, soluteMass, totalMass;
const int adjustRegion, resetTime;
extern VecR region, solInertia;
extern const real deltaT, extTemperature;
=======================================
--- /trunk/src/sysGen.c Wed Feb 9 07:48:13 2011
+++ /trunk/src/sysGen.c Tue Jun 28 02:10:55 2011
@@ -99,7 +99,7 @@
fprintf( fPtr, "done\n" ); fflush(fPtr);
}
-/* InitVels. Setting of initial COM velocities. */
+/* InitVels. Setting of initial COM velocities. Adapted from Rapaport
[p.28-29] */
void InitVels( const real extTemperature, FILE *fPtr )
{
int n;