Something went wrong on our end
Configuration.cpp 52.04 KiB
#include "Configuration.h"
#include <cassert>
#define gpuErrchk(ans) { gpuAssert((ans), __FILE__, __LINE__); }
inline void gpuAssert(cudaError_t code, const char *file, int line, bool abort=true) {
if (code != cudaSuccess) {
fprintf(stderr,"CUDA Error: %s %s %d\n", cudaGetErrorString(code), __FILE__, line);
if (abort) exit(code);
}
}
Configuration::Configuration(const char* config_file, int simNum, bool debug) :
simNum(simNum) {
// Read the parameters.
//type_d = NULL;
kTGrid_d = NULL;
//bonds_d = NULL;
//bondMap_d = NULL;
//excludes_d = NULL;
//excludeMap_d = NULL;
//angles_d = NULL;
//dihedrals_d = NULL;
setDefaults();
readParameters(config_file);
// Get the number of particles
printf("\nCounting particles specified in the ");
if (restartCoordinates.length() > 0) {
// Read them from the restart file.
printf("restart file.\n");
num = countRestart(restartCoordinates.val());
} else {
if (readPartsFromFile) readAtoms();
if (numPartsFromFile > 0) {
// Determine number of particles from input file (PDB-style)
printf("input file.\n");
num = numPartsFromFile;
} else {
// Sum up all particles in config file
printf("configuration file.\n");
//int num0 = 0;
num = 0;
for (int i = 0; i < numParts; i++) num += part[i].num;
//num = num0;
}
} // end result: variable "num" is set
// Set the number capacity
printf("\nInitial particles: %d\n", num);
if (numCap <= 0) numCap = numCapFactor*num; // max number of particles
if (numCap <= 0) numCap = 20;
// Allocate particle variables.
pos = new Vector3[num * simNum];
type = new int[num * simNum];
serial = new int[num * simNum];
posLast = new Vector3[num * simNum];
name = new String[num * simNum];
currSerial = 0;
// Now, load the coordinates
loadedCoordinates = false;
// If we have a restart file - use it
if (restartCoordinates.length() > 0) {
loadRestart(restartCoordinates.val());
printf("Loaded %d restart coordinates from `%s'.\n", num, restartCoordinates.val()); printf("Particle numbers specified in the configuration file will be ignored.\n");
loadedCoordinates = true;
} else {
// Load coordinates from a file?
if (numPartsFromFile > 0) {
loadedCoordinates = true;
for (int i = 0; i < num; i++) {
int numTokens = partsFromFile[i].tokenCount();
// Break the line down into pieces (tokens) so we can process them individually
String* tokenList = new String[numTokens];
partsFromFile[i].tokenize(tokenList);
int currType = 0;
for (int j = 0; j < numParts; j++)
if (tokenList[2] == part[j].name)
currType = j;
for (int s = 0; s < simNum; ++s)
type[i + s*num] = currType;
serial[i] = currSerial++;
pos[i] = Vector3(atof(tokenList[3].val()),
atof(tokenList[4].val()),
atof(tokenList[5].val()));
}
delete[] partsFromFile;
partsFromFile = NULL;
} else {
// Not loading coordinates from a file
populate();
if (inputCoordinates.length() > 0) {
printf("Loading coordinates from %s ... ", inputCoordinates.val());
loadedCoordinates = loadCoordinates(inputCoordinates.val());
if (loadedCoordinates)
printf("done!\n");
}
}
}
if (readBondsFromFile) readBonds();
if (readExcludesFromFile) readExcludes();
if (readAnglesFromFile) readAngles();
if (readDihedralsFromFile) readDihedrals();
kT = temperature * 0.0019872065f; // `units "k K" "kcal_mol"`
if (temperatureGridFile.length() != 0) {
printf("\nFound temperature grid file: %s\n", temperatureGridFile.val());
tGrid = new BaseGrid(temperatureGridFile.val());
printf("Loaded `%s'.\n", temperatureGridFile.val());
printf("System size %s.\n", tGrid->getExtent().toString().val());
// TODO: ask Max Belkin what this is about and how to remove hard-coded temps
float ToSo = 1.0f / (295.0f * 4.634248239f); // 1 / (To * sigma(To))
sigmaT = new BaseGrid(*tGrid);
sigmaT->shift(-122.8305f);
sigmaT->scale(0.0269167f);
sigmaT->mult(*tGrid);
sigmaT->scale(ToSo);
kTGrid = new BaseGrid(*tGrid);
float factor = 0.0019872065f; // `units "k K" "kcal_mol"`
kTGrid->scale(factor);
// char outFile[256];
// char comment[256]; sprintf(comment,"KTGrid");
// sprintf(outFile,"kTGrid.dx");
// kTGrid->write(outFile, comment);
}
printf("\nFound %d particle types.\n", numParts);
// Load the potential grids.
printf("Loading the potential grids...\n");
for (int i = 0; i < numParts; i++) {
// Decide which type of grid is given.
String map = partGridFile[i];
int len = map.length();
if (len >= 3 && map[len-3]=='.' && map[len-2]=='d' && map[len-1]=='x') {
// A dx file. Load the old-fashioned way.
part[i].pmf = new BaseGrid(map.val());
if (partGridFileScale[i] != 1.0f) part[i].pmf->scale(partGridFileScale[i]);
part[i].meanPmf = part[i].pmf->mean();
printf("Loaded dx grid `%s'.\n", map.val());
printf("System size %s.\n", part[i].pmf->getExtent().toString().val());
} else if (len >= 4 && map[len-4]=='.' && map[len-3]=='d' && map[len-2]=='e' && map[len-1]=='f') {
// A system definition file.
String rootGrid = OverlordGrid::readDefFirst(map);
OverlordGrid* over = new OverlordGrid(rootGrid.val());
int count = over->readDef(map);
printf("Loaded system def file `%s'.\n", map.val());
printf("Found %d unique grids.\n", over->getUniqueGridNum());
printf("Linked %d subgrids.\n", count);
part[i].pmf = static_cast<BaseGrid*>(over);
part[i].meanPmf = part[i].pmf->mean();
} else {
printf("WARNING: Unrecognized gridFile extension. Must be *.def or *.dx.\n");
exit(-1);
}
if (partForceXGridFile[i].length() != 0) {
part[i].forceXGrid = new BaseGrid(partForceXGridFile[i].val());
printf("Loaded `%s'.\n", partForceXGridFile[i].val());
printf("System size %s.\n", part[i].forceXGrid->getExtent().toString().val());
}
if (partForceYGridFile[i].length() != 0) {
part[i].forceYGrid = new BaseGrid(partForceYGridFile[i].val());
printf("Loaded `%s'.\n", partForceYGridFile[i].val());
printf("System size %s.\n", part[i].forceYGrid->getExtent().toString().val());
}
if (partForceZGridFile[i].length() != 0) {
part[i].forceZGrid = new BaseGrid(partForceZGridFile[i].val());
printf("Loaded `%s'.\n", partForceZGridFile[i].val());
printf("System size %s.\n", part[i].forceZGrid->getExtent().toString().val());
}
if (partDiffusionGridFile[i].length() != 0) {
part[i].diffusionGrid = new BaseGrid(partDiffusionGridFile[i].val());
printf("Loaded `%s'.\n", partDiffusionGridFile[i].val());
printf("System size %s.\n", part[i].diffusionGrid->getExtent().toString().val());
}
if (temperatureGridFile.length() != 0) {
if (partDiffusionGridFile[i].length() != 0) {
part[i].diffusionGrid->mult(*sigmaT);
} else {
part[i].diffusionGrid = new BaseGrid(*sigmaT);
part[i].diffusionGrid->scale(part[i].diffusion);
// char outFile[256];
// char comment[256]; sprintf(comment,"Diffusion for particle type %d", i);
// sprintf(outFile,"diffusion%d.dx",i);
// part[i].diffusionGrid->write(outFile, comment);
}
}
}
// Load reservoir files if any
for (int i = 0; i < numParts; i++) {
if (partReservoirFile[i].length() != 0) {
printf("\nLoading the reservoirs for %s... \n", part[i].name.val());
part[i].reservoir = new Reservoir(partReservoirFile[i].val());
int nRes = part[i].reservoir->length();
printf("\t -> %d reservoir(s) found in `%s'.\n", nRes, partReservoirFile[i].val());
}
}
// Get the system dimensions
// from the dimensions of supplied 3D potential maps
sys = part[0].pmf;
sysDim = part[0].pmf->getExtent();
// RBTODO: clean this mess up
/* // RigidBodies... */
/* if (numRigidTypes > 0) { */
/* printf("\nCounting rigid bodies specified in the configuration file.\n"); */
/* numRB = 0; */
/* // grow list of rbs */
/* for (int i = 0; i < numRigidTypes; i++) { */
/* numRB += rigidBody[i].num; */
/* std::vector<RigidBody> tmp; */
/* for (int j = 0; j < rigidBody[i].num; j++) { */
/* tmp.push_back( new RigidBody( this, rigidBody[i] ) ); */
/* } */
/* rbs.push_back(tmp); */
/* } */
// // state data
// rbPos = new Vector3[numRB * simNum];
// type = new int[numRB * simNum];
/* } */
/* printf("Initial RigidBodies: %d\n", numRB); */
// Create exclusions from the exclude rule, if it was specified in the config file
if (excludeRule != String("")) {
int oldNumExcludes = numExcludes;
Exclude* newExcludes = makeExcludes(bonds, bondMap, num, numBonds, excludeRule, numExcludes);
if (excludes == NULL) {
excludes = new Exclude[numExcludes];
} else if (numExcludes >= excludeCapacity) {
Exclude* tempExcludes = excludes;
excludes = new Exclude[numExcludes];
for (int i = 0; i < oldNumExcludes; i++)
excludes[i] = tempExcludes[i];
delete tempExcludes;
}
for (int i = oldNumExcludes; i < numExcludes; i++)
excludes[i] = newExcludes[i - oldNumExcludes];
printf("Built %d exclusions.\n",numExcludes);
// Call compareExcludeIndex with qsort to sort the excludes by BOTH ind1 AND ind2
std::sort(excludes, excludes + numExcludes, compare());
/* Each particle may have a varying number of excludes
* excludeMap is an array with one element for each particle
* which keeps track of where a particle's excludes are stored
* in the excludes array.
* excludeMap[i].x is the index in the excludes array where the ith particle's excludes begin
* excludeMap[i].y is the index in the excludes array where the ith particle's excludes end
*/
excludeMap = new int2[numPartsFromFile];
for (int i = 0; i < numPartsFromFile; i++) { excludeMap[i].x = -1;
excludeMap[i].y = -1;
}
int currPart = -1;
int lastPart = -1;
for (int i = 0; i < numExcludes; i++) {
if (excludes[i].ind1 != currPart) {
currPart = excludes[i].ind1;
excludeMap[currPart].x = i;
if (lastPart >= 0)
excludeMap[lastPart].y = i;
lastPart = currPart;
}
}
}
// Count number of particles of each type
numPartsOfType = new int[numParts];
for (int i = 0; i < numParts; ++i) {
numPartsOfType[i] = 0;
}
for (int i = 0; i < num; ++i) {
++numPartsOfType[type[i]];
}
// Some geometric stuff that should be gotten rid of.
Vector3 buffer = (sys->getCenter() + 2.0f*sys->getOrigin())/3.0f;
initialZ = buffer.z;
// Set the initial conditions.
// Do the initial conditions come from restart coordinates?
// inputCoordinates are ignored if restartCoordinates exist.
/*
if (restartCoordinates.length() > 0) {
loadRestart(restartCoordinates.val());
printf("Loaded %d restart coordinates from `%s'.\n", num, restartCoordinates.val());
printf("Particle numbers specified in the configuration file will be ignored.\n");
} else {
// Set the particle types.
// Load coordinates from a file?
if (numPartsFromFile > 0) {
for (int i = 0; i < num; i++) {
int numTokens = partsFromFile[i].tokenCount();
// Break the line down into pieces (tokens) so we can process them individually
String* tokenList = new String[numTokens];
partsFromFile[i].tokenize(tokenList);
int currType = 0;
for (int j = 0; j < numParts; j++)
if (tokenList[2] == part[j].name)
currType = j;
type[i] = currType;
serial[i] = currSerial;
currSerial++;
pos[i] = Vector3(atof(tokenList[3].val()), atof(tokenList[4].val()), atof(tokenList[5].val()));
}
if (partsFromFile != NULL) {
delete[] partsFromFile;
partsFromFile = NULL;
}
} else if (inputCoordinates.length() > 0) {
populate();
printf("Loading coordinates from %s.\n", inputCoordinates.val());
bool loaded = loadCoordinates(inputCoordinates.val());
if (loaded)
printf("Loaded initial coordinates from %s.\n", inputCoordinates.val());
}
} */
// Get the maximum particle radius.
minimumSep = 0.0f;
for (int i = 0; i < numParts; ++i)
minimumSep = std::max(minimumSep, part[i].radius);
minimumSep *= 2.5f; // Make it a little bigger.
// Default outputEnergyPeriod
if (outputEnergyPeriod < 0)
outputEnergyPeriod = 10 * outputPeriod;
// If we are running with debug ON, ask the user which force computation to use
if (debug)
getDebugForce();
printf("\n");
switchStart = cutoff - switchLen;
if (fullLongRange == 0)
printf("Cutting off the potential from %.10g to %.10g.\n", switchStart, switchStart+switchLen);
if (fullLongRange != 0)
printf("No cell decomposition created.\n");
}
Configuration::~Configuration() {
// System state
delete[] pos;
delete[] posLast;
delete[] type;
delete[] name;
// Particle parameters
delete[] part;
delete[] partGridFile;
delete[] partGridFileScale;
delete[] partForceXGridFile;
delete[] partForceYGridFile;
delete[] partForceZGridFile;
delete[] partDiffusionGridFile;
delete[] partReservoirFile;
partRigidBodyGrid.clear();
// TODO: plug memory leaks
if (partsFromFile != NULL) delete[] partsFromFile;
if (bonds != NULL) delete[] bonds;
if (bondMap != NULL) delete[] bondMap;
if (excludes != NULL) delete[] excludes;
if (excludeMap != NULL) delete[] excludeMap;
if (angles != NULL) delete[] angles;
if (dihedrals != NULL) delete[] dihedrals;
delete[] numPartsOfType;
// Table parameters
delete[] partTableFile;
delete[] partTableIndex0;
delete[] partTableIndex1;
delete[] bondTableFile;
delete[] angleTableFile;
delete[] dihedralTableFile;
//if (type_d != NULL) {
//gpuErrchk(cudaFree(type_d));
gpuErrchk(cudaFree(sys_d));
gpuErrchk(cudaFree(kTGrid_d)); gpuErrchk(cudaFree(part_d));
//gpuErrchk(cudaFree(bonds_d));
//gpuErrchk(cudaFree(bondMap_d));
//gpuErrchk(cudaFree(excludes_d));
//gpuErrchk(cudaFree(excludeMap_d));
//gpuErrchk(cudaFree(angles_d));
//gpuErrchk(cudaFree(dihedrals_d));
//}
}
void Configuration::copyToCUDA() {
printf("Copying to GPU %d\n", GPUManager::current());
BrownianParticleType **part_addr = new BrownianParticleType*[numParts];
// Copy the BaseGrid objects and their member variables/objects
gpuErrchk(cudaMalloc(&part_d, sizeof(BrownianParticleType*) * numParts));
// TODO: The above line fails when there is not enough memory. If it fails, stop.
for (int i = 0; i < numParts; i++) {
BaseGrid *pmf = NULL, *diffusionGrid = NULL;
BrownianParticleType *b = new BrownianParticleType(part[i]);
// Copy PMF
if (part[i].pmf != NULL) {
float *val = NULL;
size_t sz = sizeof(float) * part[i].pmf->getSize();
gpuErrchk(cudaMalloc(&pmf, sizeof(BaseGrid)));
gpuErrchk(cudaMalloc(&val, sz));
gpuErrchk(cudaMemcpyAsync(val, part[i].pmf->val, sz, cudaMemcpyHostToDevice));
BaseGrid *pmf_h = new BaseGrid(*part[i].pmf);
pmf_h->val = val;
gpuErrchk(cudaMemcpy(pmf, pmf_h, sizeof(BaseGrid), cudaMemcpyHostToDevice));
pmf_h->val = NULL;
}
// Copy the diffusion grid
if (part[i].diffusionGrid != NULL) {
float *val = NULL;
size_t sz = sizeof(float) * part[i].diffusionGrid->getSize();
BaseGrid *diffusionGrid_h = new BaseGrid(*part[i].diffusionGrid);
gpuErrchk(cudaMalloc(&diffusionGrid, sizeof(BaseGrid)));
gpuErrchk(cudaMalloc(&val, sz));
diffusionGrid_h->val = val;
gpuErrchk(cudaMemcpyAsync(diffusionGrid, diffusionGrid_h, sizeof(BaseGrid),
cudaMemcpyHostToDevice));
gpuErrchk(cudaMemcpy(val, part[i].diffusionGrid->val, sz, cudaMemcpyHostToDevice));
diffusionGrid_h->val = NULL;
}
b->pmf = pmf;
b->diffusionGrid = diffusionGrid;
gpuErrchk(cudaMalloc(&part_addr[i], sizeof(BrownianParticleType)));
gpuErrchk(cudaMemcpyAsync(part_addr[i], b, sizeof(BrownianParticleType),
cudaMemcpyHostToDevice));
}
// RBTODO: moved this out of preceding loop; was that correct?
gpuErrchk(cudaMemcpyAsync(part_d, part_addr, sizeof(BrownianParticleType*) * numParts,
cudaMemcpyHostToDevice));
// kTGrid_d
kTGrid_d = NULL;
if (temperatureGridFile.length() > 0) {
gpuErrchk(cudaMalloc(&kTGrid_d, sizeof(BaseGrid)));
gpuErrchk(cudaMemcpyAsync(kTGrid_d, kTGrid, sizeof(BaseGrid), cudaMemcpyHostToDevice));
}
// type_d and sys_d
gpuErrchk(cudaMalloc(&sys_d, sizeof(BaseGrid)));
gpuErrchk(cudaMemcpyAsync(sys_d, sys, sizeof(BaseGrid), cudaMemcpyHostToDevice));
/*gpuErrchk(cudaMalloc(&type_d, sizeof(int) * num * simNum)); gpuErrchk(cudaMemcpyAsync(type_d, type, sizeof(int) * num * simNum, cudaMemcpyHostToDevice));
if (numBonds > 0) {
// bonds_d
gpuErrchk(cudaMalloc(&bonds_d, sizeof(Bond) * numBonds));
gpuErrchk(cudaMemcpyAsync(bonds_d, bonds, sizeof(Bond) * numBonds, cudaMemcpyHostToDevice));
// bondMap_d
gpuErrchk(cudaMalloc(&bondMap_d, sizeof(int2) * num));
gpuErrchk(cudaMemcpyAsync(bondMap_d, bondMap, sizeof(int2) * num, cudaMemcpyHostToDevice));
}
if (numExcludes > 0) {
// excludes_d
gpuErrchk(cudaMalloc(&excludes_d, sizeof(Exclude) * numExcludes));
gpuErrchk(cudaMemcpyAsync(excludes_d, excludes, sizeof(Exclude) * numExcludes,
cudaMemcpyHostToDevice));
// excludeMap_d
gpuErrchk(cudaMalloc(&excludeMap_d, sizeof(int2) * num));
gpuErrchk(cudaMemcpyAsync(excludeMap_d, excludeMap, sizeof(int2) * num,
cudaMemcpyHostToDevice));
}
if (numAngles > 0) {
// angles_d
gpuErrchk(cudaMalloc(&angles_d, sizeof(Angle) * numAngles));
gpuErrchk(cudaMemcpyAsync(angles_d, angles, sizeof(Angle) * numAngles,
cudaMemcpyHostToDevice));
}
if (numDihedrals > 0) {
// dihedrals_d
gpuErrchk(cudaMalloc(&dihedrals_d, sizeof(Dihedral) * numDihedrals));
gpuErrchk(cudaMemcpyAsync(dihedrals_d, dihedrals,
sizeof(Dihedral) * numDihedrals,
cudaMemcpyHostToDevice));
}*/
gpuErrchk(cudaDeviceSynchronize());
}
void Configuration::setDefaults() {
// System parameters
outputName = "out";
timestep = 1e-5f;
rigidBodyGridGridPeriod = 1;
steps = 100;
seed = 0;
inputCoordinates = "";
restartCoordinates = "";
numberFluct = 0;
numberFluctPeriod = 200;
interparticleForce = 1;
tabulatedPotential = 0;
fullLongRange = 1;
// kTGridFile = ""; // Commented out for an unknown reason
temperature = 295.0f;
temperatureGridFile = "";
coulombConst = 566.440698f/92.0f;
electricField = 0.0f;
cutoff = 10.0f;
switchLen = 2.0f;
pairlistDistance = 2.0f;
imdForceScale = 1.0f;
outputPeriod = 200;
outputEnergyPeriod = -1;
outputFormat = TrajectoryWriter::formatDcd;
currentSegmentZ = -1.0f;
numCap = 0;
decompPeriod = 10; readPartsFromFile = 0;
numPartsFromFile = 0;
partsFromFile = NULL;
readBondsFromFile = false;
numBonds = 0;
bonds = NULL;
bondMap = NULL;
numTabBondFiles = 0;
readExcludesFromFile = false;
numExcludes = 0;
excludeCapacity = 256;
excludes = NULL;
excludeMap = NULL;
excludeRule = "";
readAnglesFromFile = false;
numAngles = 0;
angles = NULL;
numTabAngleFiles = 0;
readDihedralsFromFile = false;
numDihedrals = 0;
dihedrals = NULL;
numTabDihedralFiles = 0;
// Hidden parameters
// Might be parameters later
numCapFactor = 5;
}
int Configuration::readParameters(const char * config_file) {
Reader config(config_file);
printf("Read config file %s\n", config_file);
// Get the number of particles.
const int numParams = config.length();
numParts = config.countParameter("particle");
numRigidTypes = config.countParameter("rigidBody");
// Allocate the particle variables.
part = new BrownianParticleType[numParts];
partGridFile = new String[numParts];
partGridFileScale = new float[numParts];
partForceXGridFile = new String[numParts];
partForceYGridFile = new String[numParts];
partForceZGridFile = new String[numParts];
partDiffusionGridFile = new String[numParts];
partReservoirFile = new String[numParts];
partRigidBodyGrid.resize(numParts);
// Allocate the table variables.
partTableFile = new String[numParts*numParts];
partTableIndex0 = new int[numParts*numParts];
partTableIndex1 = new int[numParts*numParts];
// Allocate rigid body types
rigidBody = new RigidBodyType[numRigidTypes];
// Set a default
for (int i = 0; i < numParts; ++i) {
partGridFileScale[i] = 1.0f;
}
int btfcap = 10;
bondTableFile = new String[btfcap];
int atfcap = 10;
angleTableFile = new String[atfcap];
int dtfcap = 10;
dihedralTableFile = new String[dtfcap];
int currPart = -1;
int currTab = -1;
int currBond = -1;
int currAngle = -1;
int currDihedral = -1;
int currRB = -1;
int partClassPart = 0;
int partClassRB = 1;
int currPartClass = -1; // 0 => particle, 1 => rigidBody
for (int i = 0; i < numParams; i++) {
String param = config.getParameter(i);
String value = config.getValue(i);
// GLOBAL
if (param == String("outputName"))
outputName = value;
else if (param == String("timestep"))
timestep = (float) strtod(value.val(), NULL);
else if (param == String("rigidBodyGridGridPeriod"))
rigidBodyGridGridPeriod = atoi(value.val());
else if (param == String("steps"))
steps = atol(value.val());
else if (param == String("seed"))
seed = atoi(value.val());
else if (param == String("inputCoordinates"))
inputCoordinates = value;
else if (param == String("restartCoordinates"))
restartCoordinates = value;
else if (param == String("temperature"))
temperature = (float) strtod(value.val(),NULL);
else if (param == String("temperatureGrid"))
temperatureGridFile = value;
else if (param == String("numberFluct"))
numberFluct = atoi(value.val());
else if (param == String("numberFluctPeriod"))
numberFluctPeriod = atoi(value.val());
else if (param == String("interparticleForce"))
interparticleForce = atoi(value.val());
else if (param == String("fullLongRange") || param == String("fullElect") )
fullLongRange = atoi(value.val());
else if (param == String("coulombConst"))
coulombConst = (float) strtod(value.val(), NULL);
else if (param == String("electricField"))
electricField = (float) strtod(value.val(), NULL);
else if (param == String("cutoff"))
cutoff = (float) strtod(value.val(), NULL);
else if (param == String("switchLen"))
switchLen = (float) strtod(value.val(), NULL);
else if (param == String("pairlistDistance"))
pairlistDistance = (float) strtod(value.val(), NULL);
else if (param == String("scaleIMDForce"))
imdForceScale = (float) strtod(value.val(), NULL);
else if (param == String("outputPeriod"))
outputPeriod = atoi(value.val());
else if (param == String("outputEnergyPeriod"))
outputEnergyPeriod = atoi(value.val());
else if (param == String("outputFormat"))
outputFormat = TrajectoryWriter::getFormatCode(value);
else if (param == String("currentSegmentZ"))
currentSegmentZ = (float) strtod(value.val(), NULL);
else if (param == String("numCap"))
numCap = atoi(value.val());
else if (param == String("decompPeriod"))
decompPeriod = atoi(value.val());
// PARTICLES
else if (param == String("particle")) {
part[++currPart] = BrownianParticleType(value);
currPartClass = partClassPart;
} else if (param == String("forceXGridFile"))
partForceXGridFile[currPart] = value;
else if (param == String("forceYGridFile"))
partForceYGridFile[currPart] = value;
else if (param == String("forceZGridFile"))
partForceZGridFile[currPart] = value;
else if (param == String("diffusionGridFile"))
partDiffusionGridFile[currPart] = value;
else if (param == String("diffusion"))
part[currPart].diffusion = (float) strtod(value.val(), NULL);
else if (param == String("charge"))
part[currPart].charge = (float) strtod(value.val(), NULL);
else if (param == String("radius"))
part[currPart].radius = (float) strtod(value.val(), NULL);
else if (param == String("eps"))
part[currPart].eps = (float) strtod(value.val(), NULL);
else if (param == String("reservoirFile"))
partReservoirFile[currPart] = value;
else if (param == String("tabulatedPotential"))
tabulatedPotential = atoi(value.val());
else if (param == String("tabulatedFile"))
readTableFile(value, ++currTab);
else if (param == String("tabulatedBondFile")) {
if (numTabBondFiles >= btfcap) {
String* temp = bondTableFile;
btfcap *= 2;
bondTableFile = new String[btfcap];
for (int j = 0; j < numTabBondFiles; j++)
bondTableFile[j] = temp[j];
delete[] temp;
}
if (readBondFile(value, ++currBond))
numTabBondFiles++;
} else if (param == String("inputParticles")) {
if (readPartsFromFile) {
printf("WARNING: More than one particle file specified. Discarding new file.\n");
} else {
partFile = value;
readPartsFromFile = true;
loadedCoordinates = true;
}
} else if (param == String("inputBonds")) {
if (readBondsFromFile) {
printf("WARNING: More than one bond file specified. Discarding new bond file.\n");
} else {
bondFile = value;
readBondsFromFile = true;
}
} else if (param == String("inputExcludes")) {
if (readExcludesFromFile) {
printf("WARNING: More than one exclude file specified. Discarding new exclude file.\n");
} else {
excludeFile = value;
readExcludesFromFile = true;
}
} else if (param == String("exclude") or param == String("exclusion")) {
excludeRule = value;
} else if (param == String("inputAngles")) {
if (readAnglesFromFile) {
printf("WARNING: More than one angle file specified. Discarding new angle file.\n");
} else {
angleFile = value;
readAnglesFromFile = true;
}
} else if (param == String("tabulatedAngleFile")) {
if (numTabAngleFiles >= atfcap) {
String* temp = angleTableFile;
atfcap *= 2;
angleTableFile = new String[atfcap];
for (int j = 0; j < numTabAngleFiles; j++) angleTableFile[j] = temp[j];
delete[] temp;
}
if (readAngleFile(value, ++currAngle))
numTabAngleFiles++;
} else if (param == String("inputDihedrals")) {
if (readDihedralsFromFile) {
printf("WARNING: More than one dihedral file specified. Discarding new dihedral file.\n");
} else {
dihedralFile = value;
readDihedralsFromFile = true;
}
} else if (param == String("tabulatedDihedralFile")) {
if (numTabDihedralFiles >= dtfcap) {
String * temp = dihedralTableFile;
dtfcap *= 2;
dihedralTableFile = new String[dtfcap];
for (int j = 0; j < numTabDihedralFiles; j++)
dihedralTableFile[j] = temp[j];
delete[] temp;
}
if (readDihedralFile(value, ++currDihedral))
numTabDihedralFiles++;
} else if (param == String("gridFileScale")) {
partGridFileScale[currPart] = (float) strtod(value.val(), NULL);
} else if (param == String("rigidBodyPotential")) {
partRigidBodyGrid[currPart].push_back(value);
}
// RIGID BODY
else if (param == String("rigidBody")) {
// part[++currPart] = BrownianParticleType(value);
rigidBody[++currRB] = RigidBodyType(value, this);
currPartClass = partClassRB;
}
else if (param == String("mass"))
rigidBody[currRB].mass = (float) strtod(value.val(), NULL);
else if (param == String("inertia"))
rigidBody[currRB].inertia = stringToVector3( value );
else if (param == String("transDamping"))
rigidBody[currRB].transDamping = stringToVector3( value );
else if (param == String("rotDamping"))
rigidBody[currRB].rotDamping = stringToVector3( value );
else if (param == String("densityGrid"))
rigidBody[currRB].addDensityGrid(value);
else if (param == String("potentialGrid"))
rigidBody[currRB].addPotentialGrid(value);
else if (param == String("densityGridScale"))
rigidBody[currRB].scaleDensityGrid(value);
else if (param == String("potentialGridScale"))
rigidBody[currRB].scalePotentialGrid(value);
else if (param == String("pmfScale"))
rigidBody[currRB].scalePMF(value);
else if (param == String("position"))
rigidBody[currRB].initPos = stringToVector3( value );
else if (param == String("orientation"))
rigidBody[currRB].initRot = stringToMatrix3( value );
else if (param == String("inputRBCoordinates"))
inputRBCoordinates = value;
// COMMON
else if (param == String("num")) {
if (currPartClass == partClassPart)
part[currPart].num = atoi(value.val());
else if (currPartClass == partClassRB)
rigidBody[currRB].num = atoi(value.val());
}
else if (param == String("gridFile")) {
if (currPartClass == partClassPart)
partGridFile[currPart] = value; else if (currPartClass == partClassRB)
rigidBody[currRB].addPMF(value);
}
// UNKNOWN
else {
printf("ERROR: Unrecognized keyword `%s'.\n", param.val());
exit(1);
}
}
// extra configuration for RB types
for (int i = 0; i < numRigidTypes; i++)
rigidBody[i].setDampingCoeffs(timestep);
return numParams;
}
Vector3 Configuration::stringToVector3(String s) {
// tokenize and return
int numTokens = s.tokenCount();
if (numTokens != 3) {
printf("ERROR: could not convert input to Vector3.\n"); // TODO improve this message
exit(1);
}
String* token = new String[numTokens];
s.tokenize(token);
Vector3 v( (float) strtod(token[0], NULL),
(float) strtod(token[1], NULL),
(float) strtod(token[2], NULL) );
return v;
}
Matrix3 Configuration::stringToMatrix3(String s) {
// tokenize and return
int numTokens = s.tokenCount();
if (numTokens != 9) {
printf("ERROR: could not convert input to Matrix3.\n"); // TODO improve this message
exit(1);
}
String* token = new String[numTokens];
s.tokenize(token);
Matrix3 m( (float) strtod(token[0], NULL),
(float) strtod(token[1], NULL),
(float) strtod(token[2], NULL),
(float) strtod(token[3], NULL),
(float) strtod(token[4], NULL),
(float) strtod(token[5], NULL),
(float) strtod(token[6], NULL),
(float) strtod(token[7], NULL),
(float) strtod(token[8], NULL) );
return m;
}
void Configuration::readAtoms() {
// Open the file
FILE* inp = fopen(partFile.val(), "r");
char line[256];
// If the particle file cannot be found, exit the program
if (inp == NULL) {
printf("ERROR: Could not open `%s'.\n", partFile.val());
bool found = true;
for (int i = 0; i < numParts; i++)
if (part[i].num == 0)
found = false;
if (!found) {
printf("ERROR: Number of particles not specified in config file.\n");
exit(1);
}
printf("Using default coordinates file\n");
return;
}
// Our particle array has a starting capacity of 256
// We will expand this later if we need to.
int capacity = 256;
numPartsFromFile = 0;
partsFromFile = new String[capacity];
indices = new int[capacity];
indices[0] = 0;
// Get and process all lines of input
while (fgets(line, 256, inp) != NULL) {
// Lines in the particle file that begin with # are comments
if (line[0] == '#') continue;
String s(line);
int numTokens = s.tokenCount();
// Break the line down into pieces (tokens) so we can process them individually
String* tokenList = new String[numTokens];
s.tokenize(tokenList);
// Legitimate ATOM input lines have 6 tokens:
// ATOM | Index | Name | X-coord | Y-coord | Z-coord
// A line without exactly six tokens should be discarded.
if (numTokens != 6) {
printf("Warning: Invalid particle file line: %s\n", line);
return;
}
// Ensure that this particle's type was defined in the config file.
// If not, discard this line.
bool found;
for (int j = 0; j < numParts; j++) {
// If this particle type exists, add a new one to the list
if (part[j].name == tokenList[2]) {
found = true;
part[j].num++;
}
}
// If the particle's type does not exist according to the config file, discard it.
if (!found) {
printf("WARNING Unknown particle type %s found and discarded.\n", tokenList[2].val());
continue;
}
// If we don't have enough room in our particle array, we need to expand it.
if (numPartsFromFile >= capacity) {
// Temporary pointers to the old arrays
String* temp = partsFromFile;
int* temp2 = indices;
// Double the capacity
capacity *= 2;
// Create pointers to new arrays which are twice the size of the old ones
partsFromFile = new String[capacity];
indices = new int[capacity];
// Copy the old values into the new arrays
for (int j = 0; j < numPartsFromFile; j++) {
partsFromFile[j] = temp[j];
indices[j] = temp2[j];
}
// delete the old arrays
delete[] temp;
delete[] temp2;
}
// Make sure the index of this particle is unique. // NOTE: The particle list is sorted by index.
bool uniqueID = true;
int key = atoi(tokenList[1].val());
int mid = 0;
// If the index is greater than the last index in the list,
// this particle belongs at the end of the list. Since the
// list is kept sorted, we know this is okay.
if (numPartsFromFile == 0 || key > indices[numPartsFromFile - 1]) {
indices[numPartsFromFile] = key;
partsFromFile[numPartsFromFile++] = line;
}
// We need to do a binary search to figure out if
// the index already exists in the list.
// The assumption is that input files SHOULD have their indices sorted in
// ascending order, so we shouldn't actually use the binary search
// or the sort (which is pretty time consuming) very often.
else {
int low = 0, high = numPartsFromFile - 1;
while (low <= high) {
mid = (int)((high - low) / 2 + low);
int curr = indices[mid];
if (curr < key) {
low = mid + 1;
} else if (curr > key) {
high = mid - 1;
} else {
// For now, particles with non-unique IDs are simply not added to the array
// Other possible approaches which are not yet implemented:
// 1: Keep track of these particles and assign them new IDs after you have
// already added all of the other particles.
// 2: Get rid of ALL particles with that ID, even the ones that have already
// been added.
printf("WARNING: Non-unique ID found: %s\n", line);
uniqueID = false;
break;
}
}
if (uniqueID) {
// Add the particle to the end of the array, then sort it.
indices[numPartsFromFile] = key;
partsFromFile[numPartsFromFile++] = line;
std::sort(indices, indices + numPartsFromFile);
std::sort(partsFromFile, partsFromFile + numPartsFromFile, compare());
}
}
}
}
void Configuration::readBonds() {
// Open the file
FILE* inp = fopen(bondFile.val(), "r");
char line[256];
// If the particle file cannot be found, exit the program
if (inp == NULL) {
printf("WARNING: Could not open `%s'.\n", bondFile.val());
printf(" This simulation will not use particle bonds.\n");
return;
}
// Our particle array has a starting capacity of 256
// We will expand this later if we need to.
int capacity = 256;
numBonds = 0;
bonds = new Bond[capacity];
// Get and process all lines of input
while (fgets(line, 256, inp) != NULL) {
// Lines in the particle file that begin with # are comments
if (line[0] == '#') continue;
String s(line);
int numTokens = s.tokenCount();
// Break the line down into pieces (tokens) so we can process them individually
String* tokenList = new String[numTokens];
s.tokenize(tokenList);
// Legitimate BOND input lines have 4 tokens:
// BOND | OPERATION_FLAG | INDEX1 | INDEX2 | FILENAME
// A line without exactly five tokens should be discarded.
if (numTokens != 5) {
printf("WARNING: Invalid bond file line: %s\n", line);
continue;
}
String op = tokenList[1];
int ind1 = atoi(tokenList[2].val());
int ind2 = atoi(tokenList[3].val());
String file_name = tokenList[4];
if (ind1 == ind2) {
printf("WARNING: Invalid bond file line: %s\n", line);
continue;
}
// If we don't have enough room in our bond array, we need to expand it.
if (numBonds+1 >= capacity) { // "numBonds+1" because we are adding two bonds to array
// Temporary pointer to the old array
Bond* temp = bonds;
// Double the capacity
capacity *= 2;
// Create pointer to new array which is twice the size of the old one
bonds = new Bond[capacity];
// Copy the old values into the new array
for (int j = 0; j < numBonds; j++)
bonds[j] = temp[j];
// delete the old array
delete[] temp;
}
// Add the bond to the bond array
// We must add it twice: Once for (ind1, ind2) and once for (ind2, ind1)
// RBTODO: add ind1/2 to exclusion list here iff op == REPLACE
/*if (op == Bond::REPLACE)
if( (int)(op) == 1)
{
printf("WARNING: Bond exclusions not implemented\n");
continue;
}*/
if (ind1 < 0 || ind1 >= num || ind2 < 0 || ind2 >=num) {
printf("ERROR: Bond file line '%s' includes invalid index\n", line);
exit(1);
}
Bond* b = new Bond(op, ind1, ind2, file_name);
bonds[numBonds++] = *b;
b = new Bond(op, ind2, ind1, file_name);
bonds[numBonds++] = *b;
delete[] tokenList;
}
// Call compareBondIndex with qsort to sort the bonds by BOTH ind1 AND ind2 std::sort(bonds, bonds + numBonds, compare());
/* Each particle may have a varying number of bonds
* bondMap is an array with one element for each particle
* which keeps track of where a particle's bonds are stored
* in the bonds array.
* bondMap[i].x is the index in the bonds array where the ith particle's bonds begin
* bondMap[i].y is the index in the bonds array where the ith particle's bonds end
*/
bondMap = new int2[num];
for (int i = 0; i < num; i++) {
bondMap[i].x = -1;
bondMap[i].y = -1;
}
int currPart = -1;
int lastPart = -1;
for (int i = 0; i < numBonds; i++) {
if (bonds[i].ind1 != currPart) {
currPart = bonds[i].ind1;
bondMap[currPart].x = i;
if (lastPart >= 0) bondMap[lastPart].y = i;
lastPart = currPart;
}
}
if (bondMap[lastPart].x > 0)
bondMap[lastPart].y = numBonds;
}
void Configuration::readExcludes()
{
// Open the file
FILE* inp = fopen(excludeFile.val(), "r");
char line[256];
// If the exclusion file cannot be found, exit the program
if (inp == NULL) {
printf("WARNING: Could not open `%s'.\n", excludeFile.val());
printf("This simulation will not use exclusions.\n");
return;
}
// Our particle array has a starting capacity of 256
// We will expand this later if we need to.
excludeCapacity = 256;
numExcludes = 0;
excludes = new Exclude[excludeCapacity];
// Get and process all lines of input
while (fgets(line, 256, inp) != NULL) {
// Lines in the particle file that begin with # are comments
if (line[0] == '#') continue;
String s(line);
int numTokens = s.tokenCount();
// Break the line down into pieces (tokens) so we can process them individually
String* tokenList = new String[numTokens];
s.tokenize(tokenList);
// Legitimate EXCLUDE input lines have 3 tokens:
// BOND | INDEX1 | INDEX2
// A line without exactly three tokens should be discarded.
if (numTokens != 3) {
printf("WARNING: Invalid exclude file line: %s\n", line);
continue;
}
int ind1 = atoi(tokenList[1].val());
int ind2 = atoi(tokenList[2].val());
if (ind1 >= num || ind2 >= num)
continue;
// If we don't have enough room in our bond array, we need to expand it.
if (numExcludes >= excludeCapacity) {
// Temporary pointer to the old array
Exclude* temp = excludes;
// Double the capacity
excludeCapacity *= 2;
// Create pointer to new array which is twice the size of the old one
excludes = new Exclude[excludeCapacity];
// Copy the old values into the new array
for (int j = 0; j < numExcludes; j++)
excludes[j] = temp[j];
// delete the old array
delete[] temp;
}
// Add the bond to the exclude array
// We must add it twice: Once for (ind1, ind2) and once for (ind2, ind1)
Exclude ex(ind1, ind2);
excludes[numExcludes++] = ex;
Exclude ex2(ind2, ind1);
excludes[numExcludes++] = ex2;
delete[] tokenList;
}
// Call compareExcludeIndex with qsort to sort the excludes by BOTH ind1 AND ind2
std::sort(excludes, excludes + numExcludes, compare());
/* Each particle may have a varying number of excludes
* excludeMap is an array with one element for each particle
* which keeps track of where a particle's excludes are stored
* in the excludes array.
* excludeMap[i].x is the index in the excludes array where the ith particle's excludes begin
* excludeMap[i].y is the index in the excludes array where the ith particle's excludes end
*/
excludeMap = new int2[num];
for (int i = 0; i < num; i++) {
excludeMap[i].x = -1;
excludeMap[i].y = -1;
}
int currPart = -1;
int lastPart = -1;
for (int i = 0; i < numExcludes; i++) {
if (excludes[i].ind1 != currPart) {
currPart = excludes[i].ind1;
assert(currPart < num);
excludeMap[currPart].x = i;
if (lastPart >= 0)
excludeMap[lastPart].y = i;
lastPart = currPart;
}
}
if (excludeMap[lastPart].x > 0)
excludeMap[lastPart].y = numExcludes;
}
void Configuration::readAngles() {
FILE* inp = fopen(angleFile.val(), "r");
char line[256];
int capacity = 256;
numAngles = 0;
angles = new Angle[capacity];
// If the angle file cannot be found, exit the program
if (inp == NULL) {
printf("WARNING: Could not open `%s'.\n", angleFile.val());
printf("This simulation will not use angles.\n"); return;
}
while(fgets(line, 256, inp)) {
if (line[0] == '#') continue;
String s(line);
int numTokens = s.tokenCount();
String* tokenList = new String[numTokens];
s.tokenize(tokenList);
// Legitimate ANGLE inputs have 5 tokens
// ANGLE | INDEX1 | INDEX2 | INDEX3 | FILENAME
// Any angle input line without exactly 5 tokens should be discarded
if (numTokens != 5) {
printf("WARNING: Invalid angle input line: %s\n", line);
continue;
}
// Discard any empty line
if (tokenList == NULL)
continue;
int ind1 = atoi(tokenList[1].val());
int ind2 = atoi(tokenList[2].val());
int ind3 = atoi(tokenList[3].val());
String file_name = tokenList[4];
//printf("file_name %s\n", file_name.val());
if (ind1 >= num or ind2 >= num or ind3 >= num)
continue;
if (numAngles >= capacity) {
Angle* temp = angles;
capacity *= 2;
angles = new Angle[capacity];
for (int i = 0; i < numAngles; i++)
angles[i] = temp[i];
delete[] temp;
}
Angle a(ind1, ind2, ind3, file_name);
angles[numAngles++] = a;
delete[] tokenList;
}
std::sort(angles, angles + numAngles, compare());
// for(int i = 0; i < numAngles; i++)
// angles[i].print();
}
void Configuration::readDihedrals() {
FILE* inp = fopen(dihedralFile.val(), "r");
char line[256];
int capacity = 256;
numDihedrals = 0;
dihedrals = new Dihedral[capacity];
// If the dihedral file cannot be found, exit the program
if (inp == NULL) {
printf("WARNING: Could not open `%s'.\n", dihedralFile.val());
printf("This simulation will not use dihedrals.\n");
return;
}
while(fgets(line, 256, inp)) {
if (line[0] == '#') continue;
String s(line);
int numTokens = s.tokenCount();
String* tokenList = new String[numTokens];
s.tokenize(tokenList);
// Legitimate DIHEDRAL inputs have 6 tokens
// DIHEDRAL | INDEX1 | INDEX2 | INDEX3 | INDEX4 | FILENAME
// Any angle input line without exactly 6 tokens should be discarded
if (numTokens != 6) {
printf("WARNING: Invalid dihedral input line: %s\n", line);
continue;
}
// Discard any empty line
if (tokenList == NULL)
continue;
int ind1 = atoi(tokenList[1].val());
int ind2 = atoi(tokenList[2].val());
int ind3 = atoi(tokenList[3].val());
int ind4 = atoi(tokenList[4].val());
String file_name = tokenList[5];
//printf("file_name %s\n", file_name.val());
if (ind1 >= num or ind2 >= num
or ind3 >= num or ind4 >= num)
continue;
if (numDihedrals >= capacity) {
Dihedral* temp = dihedrals;
capacity *= 2;
dihedrals = new Dihedral[capacity];
for (int i = 0; i < numDihedrals; ++i)
dihedrals[i] = temp[i];
delete[] temp;
}
Dihedral d(ind1, ind2, ind3, ind4, file_name);
dihedrals[numDihedrals++] = d;
delete[] tokenList;
}
std::sort(dihedrals, dihedrals + numDihedrals, compare());
// for(int i = 0; i < numDihedrals; i++)
// dihedrals[i].print();
}
void Configuration::populate() {
int pn = 0;
int p = 0;
for (int i = 0; i < num; i++) {
for (int s = 0; s < simNum; ++s)
type[i + s*num] = p;
serial[i] = currSerial++;
if (++pn >= part[p].num) {
p++;
pn = 0;
}
}
}
bool Configuration::readBondFile(const String& value, int currBond) {
int numTokens = value.tokenCount();
if (numTokens != 1) {
printf("ERROR: Invalid tabulatedBondFile: %s, numTokens = %d\n", value.val(), numTokens);
return false;
}
String* tokenList = new String[numTokens];
value.tokenize(tokenList);
if (tokenList == NULL) {
printf("ERROR: Invalid tabulatedBondFile: %s; tokenList is NULL\n", value.val());
return false;
}
bondTableFile[currBond] = tokenList[0];
printf("Tabulated Bond Potential: %s\n", bondTableFile[currBond].val() );
return true;
}
bool Configuration::readAngleFile(const String& value, int currAngle) {
int numTokens = value.tokenCount();
if (numTokens != 1) {
printf("ERROR: Invalid tabulatedAngleFile: %s, numTokens = %d\n", value.val(), numTokens);
return false;
}
String* tokenList = new String[numTokens];
value.tokenize(tokenList);
if (tokenList == NULL) {
printf("ERROR: Invalid tabulatedAngleFile: %s; tokenList is NULL\n", value.val());
return false;
}
angleTableFile[currAngle] = tokenList[0];
printf("Tabulated Angle Potential: %s\n", angleTableFile[currAngle].val() );
return true;
}
bool Configuration::readDihedralFile(const String& value, int currDihedral) {
int numTokens = value.tokenCount();
if (numTokens != 1) {
printf("ERROR: Invalid tabulatedDihedralFile: %s, numTokens = %d\n", value.val(), numTokens);
return false;
}
String* tokenList = new String[numTokens];
value.tokenize(tokenList);
if (tokenList == NULL) {
printf("ERROR: Invalid tabulatedDihedralFile: %s; tokenList is NULL\n", value.val());
return false;
}
dihedralTableFile[currDihedral] = tokenList[0];
printf("Tabulated Dihedral Potential: %s\n", dihedralTableFile[currDihedral].val() );
return true;
}
void Configuration::loadRestart(const char* file_name) {
char line[STRLEN];
FILE* inp = fopen(file_name, "r");
if (inp == NULL) {
printf("GrandBrownTown:loadRestart File `%s' does not exist\n", file_name);
exit(-1);
}
int count = 0;
while (fgets(line, STRLEN, inp) != NULL) {
// Ignore comments.
int len = strlen(line);
if (line[0] == '#') continue;
if (len < 2) continue;
String s(line);
int numTokens = s.tokenCount();
if (numTokens != 4) {
printf("GrandBrownTown:loadRestart Invalid coordinate file line: %s\n", line);
fclose(inp); exit(-1);
}
String* tokenList = new String[numTokens];
s.tokenize(tokenList);
if (tokenList == NULL) {
printf("GrandBrownTown:loadRestart Invalid coordinate file line: %s\n", line);
fclose(inp);
exit(-1);
}
int typ = atoi(tokenList[0]);
float x = (float) strtod(tokenList[1],NULL);
float y = (float) strtod(tokenList[2],NULL);
float z = (float) strtod(tokenList[3],NULL);
pos[count] = Vector3(x,y,z);
type[count] = typ;
serial[count] = currSerial;
currSerial++;
if (typ < 0 || typ >= numParts) {
printf("GrandBrownTown:countRestart Invalid particle type: %d\n", typ);
fclose(inp);
exit(-1);
}
count++;
delete[] tokenList;
}
fclose(inp);
}
bool Configuration::loadCoordinates(const char* file_name) {
char line[STRLEN];
FILE* inp = fopen(file_name, "r");
if (inp == NULL) return false;
int count = 0;
while (fgets(line, STRLEN, inp) != NULL) {
// Ignore comments.
int len = strlen(line);
if (line[0] == '#') continue;
if (len < 2) continue;
String s(line);
int numTokens = s.tokenCount();
if (numTokens != 3) {
printf("ERROR: Invalid coordinate file line: %s\n", line);
fclose(inp);
return false;
}
String* tokenList = new String[numTokens];
s.tokenize(tokenList);
if (tokenList == NULL) {
printf("ERROR: Invalid coordinate file line: %s\n", line);
fclose(inp);
return false;
}
if (count >= num) {
printf("WARNING: Too many coordinates in coordinate file %s.\n", file_name);
fclose(inp);
return true;
}
float x = (float) strtod(tokenList[0],NULL);
float y = (float) strtod(tokenList[1],NULL); float z = (float) strtod(tokenList[2],NULL);
pos[count] = Vector3(x,y,z);
count++;
delete[] tokenList;
}
fclose(inp);
if (count < num) {
printf("ERROR: Too few coordinates in coordinate file.\n");
return false;
}
return true;
}
// Count the number of atoms in the restart file.
int Configuration::countRestart(const char* file_name) {
char line[STRLEN];
FILE* inp = fopen(file_name, "r");
if (inp == NULL) {
printf("ERROR: countRestart File `%s' does not exist\n", file_name);
exit(-1);
}
int count = 0;
while (fgets(line, STRLEN, inp) != NULL) {
int len = strlen(line);
// Ignore comments.
if (line[0] == '#') continue;
if (len < 2) continue;
String s(line);
int numTokens = s.tokenCount();
if (numTokens != 4) {
printf("ERROR: countRestart Invalid coordinate file line: %s\n", line);
fclose(inp);
exit(-1);
}
String* tokenList = new String[numTokens];
s.tokenize(tokenList);
if (tokenList == NULL) {
printf("ERROR: countRestart Invalid coordinate file line: %s\n", line);
fclose(inp);
exit(-1);
}
int typ = atoi(tokenList[0]);
// float x = strtod(tokenList[1],NULL);
// float y = strtod(tokenList[2],NULL);
// float z = strtod(tokenList[3],NULL);
if (typ < 0 || typ >= numParts) {
printf("ERROR: countRestart Invalid particle type: %d\n", typ);
fclose(inp);
exit(-1);
}
count++;
delete[] tokenList;
}
fclose(inp);
return count;
}
bool Configuration::readTableFile(const String& value, int currTab) {
int numTokens = value.tokenCount('@');
if (numTokens != 3) {
printf("ERROR: Invalid tabulatedFile: %s\n", value.val()); return false;
}
String* tokenList = new String[numTokens];
value.tokenize(tokenList, '@');
if (tokenList == NULL) {
printf("ERROR: Invalid tabulatedFile: %s\n", value.val());
return false;
}
partTableIndex0[currTab] = atoi(tokenList[0]);
partTableIndex1[currTab] = atoi(tokenList[1]);
partTableFile[currTab] = tokenList[2];
printf("tabulatedPotential: %d %d %s\n", partTableIndex0[currTab],
partTableIndex1[currTab], partTableFile[currTab].val() );
delete[] tokenList;
return true;
}
void Configuration::getDebugForce() {
// Allow the user to choose which force computation to use
printf("\n");
printf("(1) ComputeFull [Default] (2) ComputeSoftcoreFull\n");
printf("(3) ComputeElecFull (4) Compute (Decomposed)\n");
printf("(5) ComputeTabulated (Decomposed) (6) ComputeTabulatedFull\n");
printf("WARNING: ");
if (tabulatedPotential) {
if (fullLongRange) printf("(6) was specified by config file\n");
else printf("(5) was specified by config file\n");
} else {
if (fullLongRange != 0) printf("(%d) was specified by config file\n", fullLongRange);
else printf("(4) was specified by config file\n");
}
char buffer[256];
int choice;
while (true) {
printf("Choose a force computation (1 - 6): ");
fgets(buffer, 256, stdin);
bool good = sscanf(buffer, "%d", &choice) && (choice >= 1 && choice <= 6);
if (good)
break;
}
switch(choice) {
case 1:
tabulatedPotential = 0;
fullLongRange = 1;
break;
case 2:
tabulatedPotential = 0;
fullLongRange = 2;
break;
case 3:
tabulatedPotential = 0;
fullLongRange = 3;
break;
case 4:
tabulatedPotential = 0;
fullLongRange = 0;
break;
case 5:
tabulatedPotential = 1;
fullLongRange = 0;
break;
case 6:
tabulatedPotential = 1;
fullLongRange = 1;
break; default:
tabulatedPotential = 0;
fullLongRange = 1;
break;
}
printf("\n");
}
//////////////////////////
// Comparison operators //
//////////////////////////
bool Configuration::compare::operator()(const String& lhs, const String& rhs) {
String* list_lhs = new String[lhs.tokenCount()];
String* list_rhs = new String[rhs.tokenCount()];
lhs.tokenize(list_lhs);
rhs.tokenize(list_rhs);
int key_lhs = atoi(list_lhs[1].val());
int key_rhs = atoi(list_rhs[1].val());
delete[] list_lhs;
delete[] list_rhs;
return key_lhs < key_rhs;
}
bool Configuration::compare::operator()(const Bond& lhs, const Bond& rhs) {
int diff = lhs.ind1 - rhs.ind1;
if (diff != 0)
return lhs.ind1 < rhs.ind1;
return lhs.ind2 < rhs.ind2;
}
bool Configuration::compare::operator()(const Exclude& lhs, const Exclude& rhs) {
int diff = lhs.ind1 - rhs.ind1;
if (diff != 0)
return lhs.ind1 < rhs.ind1;
return lhs.ind2 < rhs.ind2;
}
bool Configuration::compare::operator()(const Angle& lhs, const Angle& rhs) {
int diff = lhs.ind1 - rhs.ind1;
if (diff != 0)
return lhs.ind1 < rhs.ind1;
diff = lhs.ind2 - rhs.ind2;
if (diff != 0)
return lhs.ind2 < rhs.ind2;
return lhs.ind3 < rhs.ind3;
}
bool Configuration::compare::operator()(const Dihedral& lhs, const Dihedral& rhs) {
int diff = lhs.ind1 - rhs.ind1;
if (diff != 0)
return lhs.ind1 < rhs.ind1;
diff = lhs.ind2 - rhs.ind2;
if (diff != 0)
return lhs.ind2 < rhs.ind2;
diff = lhs.ind3 - rhs.ind3;
if (diff != 0)
return lhs.ind3 < rhs.ind3;
return lhs.ind4 < rhs.ind4;
}