[brahms-md] r60 committed - More work on LAMMPS interface, introduced 'refAngle_orient' variable.
brah...@googlecode.com
Author: orsimario
Date: Mon Oct 17 09:15:55 2011
Log: More work on LAMMPS interface, introduced 'refAngle_orient'
variable.
http://code.google.com/p/brahms-md/source/detail?r=60
Modified:
/trunk/src/forceBdMain.c
/trunk/src/forceNbdMain.c
/trunk/src/in_namelist.c
/trunk/src/lammps.c
/trunk/src/main.c
/trunk/src/output.c
/trunk/src/startUp.c
/trunk/src/stepMain.c
=======================================
--- /trunk/src/forceBdMain.c Thu Aug 18 10:23:37 2011
+++ /trunk/src/forceBdMain.c Mon Oct 17 09:15:55 2011
@@ -34,6 +34,7 @@
real refAngle_PhosGlyEst; /* input parameter: reference angle for the
phosphate-glycerol-ester triplet */
real refAngle_cisUnsat; /* reference angle for tail sites corresponding to
triplets involving a cis unsaturated double bond */
real angleRigidity; /* input parameter: rigidity of the angle potential
for all angles */
+real refAngle_orient; /* input parameter: reference angle for dipole
orientation */
real orientRigidity; /* input parameter: orientational rigidity for
dipolar lipid sites */
/* CalcAngleBendInteraction -- calculates angle bending interaction,
adapted from [Rapaport, pp.283-285] */
=======================================
--- /trunk/src/forceNbdMain.c Thu Aug 18 10:23:37 2011
+++ /trunk/src/forceNbdMain.c Mon Oct 17 09:15:55 2011
@@ -63,7 +63,7 @@
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_SF( &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 );
@@ -100,8 +100,10 @@
VVAdd( siteTorque[ i ], torque );
break;
case DIP_DIP: /* dipole-dipole */
- CalcDipDipInteraction_CS( &ijPotEn, &ijForce, &ijTorque, &jiTorque,
productElectroMoments[site[i].type][site[j].type],
+ CalcDipDipInteraction_CS( &ijPotEn, &ijForce, &ijTorque, &jiTorque,
productElectroMoments[site[i].type][site[j].type],
siteOrientVec[i], siteOrientVec[j], rCutLipLip, rij, rr );
+/* CalcDipDipInteraction_SF( &ijPotEn, &ijForce, &ijTorque,
&jiTorque, productElectroMoments[site[i].type][site[j].type], */
+/* siteOrientVec[i], siteOrientVec[j], rCutLipLip, rij, rr ); */
potEnDipSum += ijPotEn;
VVAdd( ijForceTot, ijForce );
VVAdd( siteTorque[ i ], ijTorque );
=======================================
--- /trunk/src/in_namelist.c Wed Sep 28 10:00:59 2011
+++ /trunk/src/in_namelist.c Mon Oct 17 09:15:55 2011
@@ -55,7 +55,7 @@
extern real hBondWatEst, hBondWatPhos, hBondWatGly;
extern real hBondAmineWat, hBondAminePhos, hBondAmineEst;
extern real epsTail, epsEst, sigTail, angleRigidity, orientRigidity;
-extern real refAngle_cisUnsat, refAngle_PhosGlyEst, refAngle_CholPhosGly;
+extern real refAngle_cisUnsat, refAngle_PhosGlyEst, refAngle_CholPhosGly,
refAngle_orient;
extern real massChol, massAmine, massPhos, massGly, massEst, massTail,
massWat;
extern real inertiaWat, inertiaEst, inertiaGly;
extern real cholCharge, phosCharge, glyDipole, estDipole, watDipole,
rangeWatRdf;
@@ -85,6 +85,7 @@
refAngle_PhosGlyEst *= degrees_IN_radians; /* convert to radians */
refAngle_CholPhosGly *= degrees_IN_radians; /* convert to radians */
refAngle_cisUnsat *= degrees_IN_radians; /* convert to radians */
+ refAngle_orient *= degrees_IN_radians; /* convert to radians */
fprintf(fp, "done\n"); fflush(fp);
}
@@ -163,10 +164,11 @@
NameR( estDipole ),
NameR( springRigidity ),
NameR( scaleRefSpringLength ),
- NameR( angleRigidity ),
NameR( refAngle_CholPhosGly ),
NameR( refAngle_PhosGlyEst ),
NameR( refAngle_cisUnsat ),
+ NameR( angleRigidity ),
+ NameR( refAngle_orient ),
NameR( orientRigidity ),
NameR( massChol ),
NameR( massAmine ),
=======================================
--- /trunk/src/lammps.c Wed Sep 28 10:00:59 2011
+++ /trunk/src/lammps.c Mon Oct 17 09:15:55 2011
@@ -14,9 +14,10 @@
extern const VecR region, *siteCoords, *siteOrientVec;
extern const Site *site;
-extern const real timeNow, *mass, *charge, *dipole, *sigma, **sig, **eps,
massAmine;
+extern const real rCutLipLip;
+extern const real timeNow, *mass, *charge, *dipole, *sigma, **sig, **eps,
massAmine, scaleRefSpringLength, refAngle_orient;
extern const real springRigidity, angleRigidity, orientRigidity,
refAngle_CholPhosGly, refAngle_PhosGlyEst, refAngle_cisUnsat;
-extern const int nSites, nLipids, nDOPCsDSPCs, nDOPEs, nTypes, nWaters;
+extern const int nSites, nLipids, nDOPCsDSPCs, nDOPEs, nTypes, nWaters,
excludeBondedNebrs;
extern int zReplicas;
// local variables:
@@ -41,7 +42,6 @@
static const int _cisUnsat_ = 4;
static const int DIPOLE = 5;
-
/*******************************************************************************************************
* WriteLammpsData -- writes "data." file to be read into the LAMMPS
program through read_data command *
*******************************************************************************************************/
@@ -222,6 +222,37 @@
}
if (angleCount!=nAngles) {printf("Error:
angleCount=%d!=%d=nAngles.\n",angleCount,nAngles); exit(0);}
}
+
+/*******************************************************************************************************
+ * DumpLammpsData -- dumps "data." file to be read into the LAMMPS program
through read_data command *
+
*******************************************************************************************************/
+void DumpLammpsData()
+{
+ FILE *fp; // file pointer
+ char fn[20]; // file name
+ char nString[10];
+
+ fn[ 0 ] = '\0';
+ strcat( fn, "data." );
+
+ if (nLipids) {
+ // if (nDOPCsDSPCs) {
+ sprintf( nString, "%d", nLipids );
+ strcat( fn, nString );
+ strcat( fn, "lip");
+ // }
+ }
+
+ if (nWaters) {
+ sprintf( nString, "%d", nWaters );
+ strcat( fn, nString );
+ strcat( fn, "wat");
+ }
+
+ fp = fopen( fn, "w" );
+ WriteLammpsData( fp );
+ fclose(fp);
+}
/**********************************************************************************************
* WriteLammpsForcefield -- writes "forcefield" file to include into
LAMMPS input script *
@@ -235,7 +266,12 @@
real radians_in_degrees = ( 180 / PI ); // conversion factor
// fprintf(fp, "pair_style hybrid lj/cut 12 lj/cut/coul/cut 12
dipole/cut 9\n\n");
- fprintf(fp, "pair_style dipole/sf 12\n");
+ fprintf(fp, "pair_style dipole/sf %2.0f\n", rCutLipLip*LENGTH_Angstrom);
+
+ if (excludeBondedNebrs)
+ fprintf(fp, "special_bonds lj/coul 0.0 1.0 1.0\n");
+ else fprintf(fp, "special_bonds lj/coul 1.0 1.0 1.0\n");
+
fprintf(fp, "bond_style harmonic\n");
fprintf(fp, "angle_style hybrid cosine/squared dipole\n");
fprintf(fp, "\n");
@@ -287,20 +323,20 @@
// syntax: bond_coeff bondType rigidityConstant referenceLength
// NOTE: in LAMMPS, this potential is defined as K(r-r0)^2, so
K_LAMMPS = k_brahms/2
if (nLipids==nDOPCsDSPCs)
- fprintf(fp, "bond_coeff%5d %7.3f%7.3f\n", X_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", X_PHOS,
springRigidity_LAMMPS_units/2,
scaleRefSpringLength*sig[CHOLINE_TYPE][PHOSPHATE_TYPE]*nm_IN_Angstrom);
else if (nLipids==nDOPEs)
- fprintf(fp, "bond_coeff%5d %7.3f%7.3f\n", X_PHOS,
springRigidity_LAMMPS_units/2,
0.9*sig[AMINE_TYPE][PHOSPHATE_TYPE]*nm_IN_Angstrom);
+ fprintf(fp, "bond_coeff%5d %7.3f%7.3f\n", X_PHOS,
springRigidity_LAMMPS_units/2,
scaleRefSpringLength*sig[AMINE_TYPE][PHOSPHATE_TYPE]*nm_IN_Angstrom);
else {printf("Error: bond assignment for mixed lipid bilayers to be
coded.\n"); exit(0);}
- 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", 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, "bond_coeff%5d %7.3f%7.3f\n", PHOS_GLY,
springRigidity_LAMMPS_units/2,
scaleRefSpringLength*sig[PHOSPHATE_TYPE][GLYCEROL_TYPE]*nm_IN_Angstrom);
+ fprintf(fp, "bond_coeff%5d %7.3f%7.3f\n", GLY_EST,
springRigidity_LAMMPS_units/2,
scaleRefSpringLength*sig[GLYCEROL_TYPE][ESTER_TYPE]*nm_IN_Angstrom);
+ fprintf(fp, "bond_coeff%5d %7.3f%7.3f\n", EST_CHAIN,
springRigidity_LAMMPS_units/2,
scaleRefSpringLength*sig[ESTER_TYPE][TAIL_TYPE]*nm_IN_Angstrom);
+ fprintf(fp, "bond_coeff%5d %7.3f%7.3f\n", CHAIN_CHAIN,
springRigidity_LAMMPS_units/2,
scaleRefSpringLength*sig[TAIL_TYPE][TAIL_TYPE]*nm_IN_Angstrom);
+ // fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", AMINE_PHOS,
springRigidity_LAMMPS_units/2,
scaleRefSpringLength*sig[PHOSPHATE_TYPE][AMINE_TYPE]*nm_IN_Angstrom);
+ /* fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 10,
springRigidity_LAMMPS_units/2,
scaleRefSpringLength*sig[][]*nm_IN_Angstrom); */
+ /* fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 11,
springRigidity_LAMMPS_units/2,
scaleRefSpringLength*sig[][]*nm_IN_Angstrom); */
+ /* fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 12,
springRigidity_LAMMPS_units/2,
scaleRefSpringLength*sig[][]*nm_IN_Angstrom); */
+ /* fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 13,
springRigidity_LAMMPS_units/2,
scaleRefSpringLength*sig[][]*nm_IN_Angstrom); */
+ /* fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 14,
springRigidity_LAMMPS_units/2,
scaleRefSpringLength*sig[][]*nm_IN_Angstrom); */
fprintf(fp, "\n");
// angle_coeff for 'cosine/squared' style
@@ -310,7 +346,23 @@
fprintf(fp, "angle_coeff%5d cosine/squared%8.3f%8.2f\n", _gly_,
angleRigidity_LAMMPS_units/2, refAngle_PhosGlyEst*radians_in_degrees);
fprintf(fp, "angle_coeff%5d cosine/squared%8.3f%8.2f\n", _piSat_,
angleRigidity_LAMMPS_units/2, PI*radians_in_degrees);
fprintf(fp, "angle_coeff%5d cosine/squared%8.3f%8.2f\n", _cisUnsat_,
angleRigidity_LAMMPS_units/2, refAngle_cisUnsat*radians_in_degrees);
- fprintf(fp, "angle_coeff%5d dipole %8.3f%8.2f\n", DIPOLE,
orientRigidity_LAMMPS_units/2, 0.); // '0' is irrelevant at the moment, as
gamma_reference = 0
+ fprintf(fp, "angle_coeff%5d dipole %8.3f%8.2f\n", DIPOLE,
orientRigidity_LAMMPS_units/2, refAngle_orient*radians_in_degrees);
fprintf(fp, "\n");
}
}
+
+/*******************************************************************************************
+ * DumpLammpsForcefield -- dumps "forcefield." file to be included into
the LAMMPS program *
+
*******************************************************************************************/
+void DumpLammpsForcefield()
+{
+ FILE *fp; // file pointer
+ char fn[20]; // file name
+
+ fn[ 0 ] = '\0';
+ strcat( fn, "forcefield.elba" );
+
+ fp = fopen( fn, "w" );
+ WriteLammpsForcefield( fp );
+ fclose(fp);
+}
=======================================
--- /trunk/src/main.c Wed Sep 28 10:00:59 2011
+++ /trunk/src/main.c Mon Oct 17 09:15:55 2011
@@ -201,6 +201,8 @@
void ComputeAndUnwrapSoluteLatCom();
void ComputeLipLatCom();
void ConvertInputParametersToBrahmsUnits();
+void DumpLammpsData();
+void DumpLammpsForcefield();
void EvalDiffusion();
void EvalEdp();
void EvalLatDiff();
@@ -279,13 +281,8 @@
SetupJob( fPtr );
fclose( fPtr );
- fPtr = fopen( "forcefield.elba", "w" );
- WriteLammpsForcefield( fPtr );
- fclose( fPtr );
-
- fPtr = fopen( "data.initial", "w" );
- WriteLammpsData( fPtr );
- fclose( fPtr );
+ DumpLammpsForcefield( fPtr );
+ DumpLammpsData( fPtr );
SetupInterrupt();
moreCycles = 1;
=======================================
--- /trunk/src/output.c Wed Sep 28 10:00:59 2011
+++ /trunk/src/output.c Mon Oct 17 09:15:55 2011
@@ -34,7 +34,7 @@
void PrintBrahmsVersion()
{
printf("\n**********************\n");
- printf( "* BRAHMS - 27Sep2011 *\n");
+ printf( "* BRAHMS - 17Oct2011 *\n");
printf( "**********************\n\n");
fflush(stdout);
}
=======================================
--- /trunk/src/startUp.c Wed Sep 28 10:00:59 2011
+++ /trunk/src/startUp.c Mon Oct 17 09:15:55 2011
@@ -474,191 +474,6 @@
printf( "then every %.1f ns.\n", stepAvgLatDiff * deltaT * TIME_ns );
InitLatDiff();
}
-
-static void InsertSoluteForZConstraint()
-{
- FILE *fPtr;
- // the solute is the last entry in the site array, i.e., site[ nSites -
1 ]
- if ( 0 == nSolutes ) { // special case: the "solute" is actually an SSD
site!
- site[ nSites - 1 ].mechType = NSRB;
- site[ nSites - 1 ].type = WATER_TYPE;
- site[ nSites - 1 ].lipidUnit = WATER;
- } else { // general case
- site[ nSites - 1 ].mechType = NSRB;
- site[ nSites - 1 ].type = SOLUTE_TYPE;
- site[ nSites - 1 ].lipidUnit = SOLUTE;
- }
- // set the rotation matrix equal to another water, e.g., nsites-2
- MCopy( siteRotMat[ nSites - 1 ].u, siteRotMat[ nSites - 2 ].u );
- printf("*** Solute coords before setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-1].x, siteCoords[nSites-1].y, siteCoords[nSites-1].z);
-
- if ( ( fPtr = fopen( "zConstraintDistanceSolute.in", "r" ) ) == 0 ) {
- printf( "*** error: Can't open file *zConstraintDistanceSolute.in* -
make sure it's here!\nExiting...\n" );
- exit( 0 );
- } else {
- fscanf( fPtr, "%lf", &zConstraintDist );
- if (zConstraintDist>0) printf("*** Warning!!! zConstraintDist is
positive - normally it should be set as negative!\n");
- }
-
- VSet( siteCoords[ nSites - 1 ], 0., 0., zConstraintDist/LENGTH_nm ); //
setting and scaling
- printf("*** Solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-1].x, siteCoords[nSites-1].y, siteCoords[nSites-1].z);
- VZero( siteVelocity[ nSites - 1 ] ); // initialize solute's velocity
- VZero( siteForce[ nSites - 1 ] ); // initialize force on solute
- VZero( siteAngMom[ nSites - 1 ] ); // initialize solute's angular
momentum
- VZero( siteTorque[ nSites - 1 ] ); // initialize torque on solute
-}
-
-static void InsertMultipleSolutes()
-{
- // the solute is the last entry in the site array, i.e., site[ nSites -
1 ]
- // for 27/64 solutes ref: rapaport p.28-31) */
-
- int i, ix, iy, iz;
- real d;
- VecI uCell;
- VecR gap, c;
-
- for ( i = 1; i <= nSolutes; i++ ) {
-
- site[ nSites - i ].mechType = NSRB;
- site[ nSites - i ].type = SOLUTE_TYPE;
- site[ nSites - i ].lipidUnit = SOLUTE;
-
- // set the rotation matrix equal to another water, e.g., nsites - (
nSolutes + i )
- MCopy( siteRotMat[ nSites - i ].u, siteRotMat[ nSites - ( nSolutes + i
) ].u );
- }
-
- VSet( uCell, 0, 0, 0 );
-
- switch ( nSolutes ) {
- case 27:
- VSet( uCell, 3, 3, 3 );
- break;
- case 45:
- VSet( uCell, 3, 3, 5 );
- break;
- case 64:
- VSet( uCell, 4, 4, 4 );
- break;
- }
-
- // set solutes coordinates
- switch ( nSolutes ) {
- case 1:
- VSet( siteCoords[nSites-1], 0, 0, 0 );
- printf("*** solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-1].x, siteCoords[nSites-1].y, siteCoords[nSites-1].z);
- break;
- case 2:
- d = 5;
- // VSet( siteCoords[nSites-1], region.x / 4, region.y / 4,
region.z / 4 );
- VSet( siteCoords[nSites-1], d, d, d );
- printf("*** solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-1].x, siteCoords[nSites-1].y, siteCoords[nSites-1].z);
- // VSet( siteCoords[nSites-2],-region.x / 4,-region.y / 4,-region.z
/ 4 );
- VSet( siteCoords[nSites-2],-d,-d,-d );
- printf("*** solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-2].x, siteCoords[nSites-2].y, siteCoords[nSites-2].z);
- break;
- case 4:
- d = 5;
- // VSet( siteCoords[nSites-1], region.x / 4, region.y / 4,
region.z / 4 );
- VSet( siteCoords[nSites-1], d, d, d );
- printf("*** solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-1].x, siteCoords[nSites-1].y, siteCoords[nSites-1].z);
- // VSet( siteCoords[nSites-2],-region.x / 4,-region.y / 4,-region.z
/ 4 );
- VSet( siteCoords[nSites-2],-d,-d,-d );
- printf("*** solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-2].x, siteCoords[nSites-2].y, siteCoords[nSites-2].z);
- VSet( siteCoords[nSites-3], d,-d, d );
- printf("*** solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-3].x, siteCoords[nSites-3].y, siteCoords[nSites-3].z);
- VSet( siteCoords[nSites-4],-d, d,-d );
- printf("*** solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-4].x, siteCoords[nSites-4].y, siteCoords[nSites-4].z);
- break;
- case 6:
- d = 5;
- VSet( siteCoords[nSites-1], d, d, d );
- printf("*** solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-1].x, siteCoords[nSites-1].y, siteCoords[nSites-1].z);
- VSet( siteCoords[nSites-2],-d,-d,-d );
- printf("*** solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-2].x, siteCoords[nSites-2].y, siteCoords[nSites-2].z);
- VSet( siteCoords[nSites-3], d,-d, d );
- printf("*** solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-3].x, siteCoords[nSites-3].y, siteCoords[nSites-3].z);
- VSet( siteCoords[nSites-4],-d, d,-d );
- printf("*** solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-4].x, siteCoords[nSites-4].y, siteCoords[nSites-4].z);
- VSet( siteCoords[nSites-5], d,-d,-d );
- printf("*** solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-5].x, siteCoords[nSites-5].y, siteCoords[nSites-5].z);
- VSet( siteCoords[nSites-6],-d,-d, d );
- break;
- case 8:
- d = 5;
- VSet( siteCoords[nSites-1], d, d, d );
- printf("*** solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-1].x, siteCoords[nSites-1].y, siteCoords[nSites-1].z);
- VSet( siteCoords[nSites-2],-d,-d,-d );
- printf("*** solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-2].x, siteCoords[nSites-2].y, siteCoords[nSites-2].z);
- VSet( siteCoords[nSites-3], d,-d, d );
- printf("*** solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-3].x, siteCoords[nSites-3].y, siteCoords[nSites-3].z);
- VSet( siteCoords[nSites-4],-d, d,-d );
- printf("*** solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-4].x, siteCoords[nSites-4].y, siteCoords[nSites-4].z);
- VSet( siteCoords[nSites-5], d,-d,-d );
- printf("*** solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-5].x, siteCoords[nSites-5].y, siteCoords[nSites-5].z);
- VSet( siteCoords[nSites-6],-d,-d, d );
- printf("*** solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-6].x, siteCoords[nSites-6].y, siteCoords[nSites-6].z);
- VSet( siteCoords[nSites-7], d, d,-d );
- printf("*** solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-7].x, siteCoords[nSites-7].y, siteCoords[nSites-7].z);
- VSet( siteCoords[nSites-8],-d, d, d );
- printf("*** solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-8].x, siteCoords[nSites-8].y, siteCoords[nSites-8].z);
- break;
- case 16:
- d = 5;
- VSet( siteCoords[nSites-1], d, d, d );
- printf("*** solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-1].x, siteCoords[nSites-1].y, siteCoords[nSites-1].z);
- VSet( siteCoords[nSites-2],-d,-d,-d );
- printf("*** solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-2].x, siteCoords[nSites-2].y, siteCoords[nSites-2].z);
- VSet( siteCoords[nSites-3], d,-d, d );
- printf("*** solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-3].x, siteCoords[nSites-3].y, siteCoords[nSites-3].z);
- VSet( siteCoords[nSites-4],-d, d,-d );
- printf("*** solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-4].x, siteCoords[nSites-4].y, siteCoords[nSites-4].z);
- VSet( siteCoords[nSites-5], d,-d,-d );
- printf("*** solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-5].x, siteCoords[nSites-5].y, siteCoords[nSites-5].z);
- VSet( siteCoords[nSites-6],-d,-d, d );
- printf("*** solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-6].x, siteCoords[nSites-6].y, siteCoords[nSites-6].z);
- VSet( siteCoords[nSites-7], d, d,-d );
- printf("*** solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-7].x, siteCoords[nSites-7].y, siteCoords[nSites-7].z);
- VSet( siteCoords[nSites-8],-d, d, d );
- printf("*** solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-8].x, siteCoords[nSites-8].y, siteCoords[nSites-8].z);
- VSet( siteCoords[nSites-9], d/2., d/2., d/2. );
- printf("*** solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-9].x, siteCoords[nSites-9].y, siteCoords[nSites-9].z);
- VSet( siteCoords[nSites-10],-d/2.,-d/2.,-d/2. );
- printf("*** solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-10].x, siteCoords[nSites-10].y, siteCoords[nSites-10].z);
- VSet( siteCoords[nSites-11], d/2.,-d/2., d/2. );
- printf("*** solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-11].x, siteCoords[nSites-11].y, siteCoords[nSites-11].z);
- VSet( siteCoords[nSites-12],-d/2., d/2.,-d/2. );
- printf("*** solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-12].x, siteCoords[nSites-12].y, siteCoords[nSites-12].z);
- VSet( siteCoords[nSites-13], d/2.,-d/2.,-d/2. );
- printf("*** solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-13].x, siteCoords[nSites-13].y, siteCoords[nSites-13].z);
- VSet( siteCoords[nSites-14],-d/2.,-d/2., d/2. );
- printf("*** solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-14].x, siteCoords[nSites-14].y, siteCoords[nSites-14].z);
- VSet( siteCoords[nSites-15], d/2., d/2.,-d/2. );
- printf("*** solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-15].x, siteCoords[nSites-15].y, siteCoords[nSites-15].z);
- VSet( siteCoords[nSites-16],-d/2., d/2., d/2. );
- printf("*** solute coords after setting [nm]: %.3f, %.3f, %.3f\n",
siteCoords[nSites-16].x, siteCoords[nSites-16].y, siteCoords[nSites-16].z);
- break;
- case 64: // fall through
- case 45: // fall through
- case 27:
- VDiv( gap, region, uCell ); // gap is the edge of each subcell that
will contain a solute
- i = 1;
- for ( iz = 0; iz < uCell.z; iz++ ) {
- for ( iy = 0; iy < uCell.y; iy++ ) {
- for ( ix = 0; ix < uCell.x; ix++ ) {
- VSet( c, ix + .5, iy + .5, iz + .5 );
- VMul ( c, c, gap );
- VVSAdd( c, -.5, region );
- siteCoords[ nSites - i ] = c;
- printf( "Assigning coordinates to solute # %d\n", i );
- ++i;
- }
- }
- }
- break;
- default: printf( "ERROR in InsertSolute(): the number of solutes can
only be 2,4,6, 8, 16, 27, 45 or 64! Exiting...\n" ); exit( 0 );
- }
-}
static void FindSoluteNAtoms() // open solute.in file and count lines to
work out the number of atoms
{
@@ -727,24 +542,6 @@
totalMass += mass[ site[ n ].type ];
}
}
-
-static void ConvertSoluteUnits()
-{
- int a;
-
- printf("### Converting solute units - reminder: review this before
reimplementing AL-CG functionalities!!! Exiting...\n");
- exit(0);
-
- soluteMass /= MASS_amu;
- VScale( solInertia, 1 / ( Sqr( LENGTH_Angstrom ) * MASS_amu ) );
-
- DO_ATOM {
- VScale( atomCoordinates[ a ], 1. / LENGTH_Angstrom );
- atom[ a ].charge *= e_IN_BRAHMS_UNITS;
- atom[ a ].sig /= LENGTH_Angstrom;
- atom[ a ].eps /= BRAHMS_ENERGY_IN_kcal_PER_mol;
- }
-}
void SetupJob( FILE *fp) // Adapted from [Rapaport, pp.20-30]
{
=======================================
--- /trunk/src/stepMain.c Wed Sep 28 10:00:59 2011
+++ /trunk/src/stepMain.c Mon Oct 17 09:15: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.
@@ -257,78 +257,6 @@
//} else VMul( siteCoords[ n ], siteCoords[ n ], region );
}
}
-
-static void LoadSoluteParametersFromFile()
-{
- FILE *fPtr;
- char line[ 81 ], *token;
- double d;
- int a;
-
- if ( ( fPtr = fopen( "solute.in", "r" ) ) == 0 ) {
- printf( "*** Error: Can't open file *solute.in* - check that it is
actually here!\nExiting...\n" );
- exit( 0 );
- }
- DO_ATOM {
- fgets( line, 81, fPtr );
- token = strtok( line, " \t\n" ); // find first token
- strcpy( atom[ a ].type, token );
- token = strtok( NULL, " \t\n" ); //
- d = atof( token );
- atomCoordinates[ a ].x = d;
- token = strtok( NULL, " \t\n" ); //
- d = atof( token );
- atomCoordinates[ a ].y = d;
- token = strtok( NULL, " \t\n" ); //
- d = atof( token );
- atomCoordinates[ a ].z = d;
- token = strtok( NULL, " \t\n" ); //
- d = atof( token );
- atom[ a ].charge = d;
- token = strtok( NULL, " \t\n" ); //
- d = atof( token );
- atom[ a ].sig = d;
- token = strtok( NULL, " \t\n" ); //
- d = atof( token );
- atom[ a ].eps = d;
- }
- fclose( fPtr );
-}
-
-static void ConvertSoluteParameters()
-{
- int a;
- DO_ATOM {
- VScale( atomCoordinates[ a ], Angstrom_IN_BRAHMS_UNITS );
- atom[ a ].charge *= e_IN_BRAHMS_UNITS;
- atom[ a ].sig *= Angstrom_IN_BRAHMS_UNITS;
- atom[ a ].eps /= BRAHMS_ENERGY_IN_kcal_PER_mol;
- }
-}
-
-static void ScaleUpSoluteParameters( real soluteGrowingFactor )
-{
- int a;
- DO_ATOM {
- VScale( atomCoordinates[ a ], soluteGrowingFactor );
- atom[ a ].charge *= soluteGrowingFactor;
- atom[ a ].sig *= soluteGrowingFactor;
- atom[ a ].eps *= soluteGrowingFactor;
- }
-}
-
-static void GrowSolute()
-{ // slowly scales up the solute(s)' size and parameters from 0 to the
actual values, step by step (as it's called at every step)
- real solGrowFac = ( real ) stepCount / stepLimit; // growing factor,
slowly going from 0 to 1 over the growing run - remember to set stepLimit
properly!
- //printf("*** solGrowFac = %f\n", solGrowFac);
- // printf("*** Step %4d LoadSoluteParametersFromFile();\t", stepCount);
- LoadSoluteParametersFromFile();
- // printf("*** ConvertSoluteParameters();\t", stepCount);
- ConvertSoluteParameters();
- // printf("*** ScaleUpSoluteParameters... \t", stepCount);
- ScaleUpSoluteParameters( solGrowFac );
- // printf("*** ScaleUpSoluteParameters done; \n", stepCount);
-}
/********************************************************************************
* StoreWatCoordsPrevStep() -- store water coordinates at the previous
timestep *