From 6f782b6d68d1622657df98d99f5656cf92119d1f Mon Sep 17 00:00:00 2001
From: Chris Maffeo <cmaffeo2@illinois.edu>
Date: Tue, 11 Aug 2020 13:27:18 -0500
Subject: [PATCH] Added 4-particle vecangle potential; untested, and
configuration.cpp still needs to build the potential dictionary using
type,filename tuples as keys
---
src/Configuration.cpp | 45 +++++++++++++------
src/TabulatedPotential.h | 93 +++++++++++++++++++++++++++++-----------
2 files changed, 100 insertions(+), 38 deletions(-)
diff --git a/src/Configuration.cpp b/src/Configuration.cpp
index b08fdd3..5b3412a 100644
--- a/src/Configuration.cpp
+++ b/src/Configuration.cpp
@@ -2048,7 +2048,7 @@ void Configuration::readProductPotentials() {
printf("\rDEBUG: reading line %d",numProductPotentials+1);
// Legitimate ProductPotential inputs have at least 7 tokens
- // BONDANGLE | INDEX1 | INDEX2 | INDEX3 | POT_FILENAME1 | INDEX4 | INDEX5 | POT_FILENAME2 ...
+ // BONDANGLE | INDEX1 | INDEX2 | INDEX3 | [TYPE1] | POT_FILENAME1 | INDEX4 | INDEX5 | [TYPE2] POT_FILENAME2 ...
if (numTokens < 7) {
printf("WARNING: Invalid product potential input line (too few tokens %d): %s\n", numTokens, line);
continue;
@@ -2058,27 +2058,44 @@ void Configuration::readProductPotentials() {
if (tokenList == NULL)
continue;
+ SimplePotentialType type = BOND; // initialize to suppress warning
+ bool type_specified = false;
for (int i = 1; i < numTokens; ++i) {
char *end;
// printf("DEBUG: Working on token %d '%s'\n", i, tokenList[i].val());
+
+ // Try to convert token to integer
int index = (int) strtol(tokenList[i].val(), &end, 10);
if (tokenList[i].val() == end || *end != '\0' || errno == ERANGE) {
- indices.push_back(tmp);
- String& n = tokenList[i];
- pot_names.push_back( n );
- if ( simple_potential_ids.find(n) == simple_potential_ids.end() ) {
- // Could not find fileName in dictionary, so read and add it
- unsigned int s = tmp.size();
- if (s < 2 || s > 4) {
- printf("WARNING: Invalid product potential input line (indices of potential %d == %d): %s\n", i, s, line);
- continue;
+ // Failed to convert token to integer; therefore it must be a potential name or type
+
+ // Try to match a type
+ String n = tokenList[i];
+ n.lower();
+ if (n == "bond") { type = DIHEDRAL; type_specified = true; }
+ else if (n == "angle") { type = DIHEDRAL; type_specified = true; }
+ else if (n == "dihedral") { type = DIHEDRAL; type_specified = true; }
+ else if (n == "vecangle") { type = VECANGLE; type_specified = true; }
+ else { // Not a type, therefore a path to a potential
+ n = tokenList[i];
+ indices.push_back(tmp);
+ pot_names.push_back( n );
+ // TODO: Key should be tuple of (type,n)
+ if ( simple_potential_ids.find(n) == simple_potential_ids.end() ) {
+ // Could not find fileName in dictionary, so read and add it
+ unsigned int s = tmp.size();
+ if (s < 2 || s > 4) {
+ printf("WARNING: Invalid product potential input line (indices of potential %d == %d): %s\n", i, s, line);
+ continue;
+ }
+ simple_potential_ids[n] = simple_potentials.size();
+ if (not type_specified) type = s==2? BOND: s==3? ANGLE: DIHEDRAL;
+ simple_potentials.push_back( SimplePotential(n.val(), type) );
}
+ tmp.clear();
+ type_specified = false;
- simple_potential_ids[n] = simple_potentials.size();
- SimplePotentialType t = s==2? BOND: s==3? ANGLE: DIHEDRAL;
- simple_potentials.push_back( SimplePotential(n.val(), t) );
}
- tmp.clear();
} else {
if (index >= num) {
continue;
diff --git a/src/TabulatedPotential.h b/src/TabulatedPotential.h
index 8c6f812..0aa73a0 100644
--- a/src/TabulatedPotential.h
+++ b/src/TabulatedPotential.h
@@ -207,7 +207,7 @@ private:
/* } */
-enum SimplePotentialType { UNSET, BOND, ANGLE, DIHEDRAL };
+enum SimplePotentialType { UNSET, BOND, ANGLE, DIHEDRAL, VECANGLE };
// enum PotentialTypeAtoms { bond=2, angle=3, dihedral=4 };
@@ -250,6 +250,8 @@ public:
val = compute_angle(pos, sys, particles[0], particles[1], particles[2]);
else if (type == DIHEDRAL)
val = compute_dihedral(pos, sys, particles[0], particles[1], particles[2], particles[3]);
+ else if (type == VECANGLE)
+ val = compute_vecangle(pos, sys, particles[0], particles[1], particles[2], particles[3]);
return val;
}
@@ -275,18 +277,25 @@ public:
const Vector3 ab = sys->wrapDiff(pos[j] - pos[i]);
const Vector3 bc = sys->wrapDiff(pos[k] - pos[j]);
const Vector3 ac = sys->wrapDiff(pos[k] - pos[i]);
+ return compute_angle( ab.length2(), bc.length2(), ac.length2() );
+ }
- // Find the distance between each pair of particles
- float distab = ab.length2();
- float distbc = bc.length2();
- const float distac2 = ac.length2();
+ HOST DEVICE inline float compute_vecangle(const Vector3* __restrict__ pos,
+ const BaseGrid* __restrict__ sys,
+ int i, int j, int k, int l) const {
+ const Vector3 ab = sys->wrapDiff(pos[j] - pos[i]);
+ const Vector3 bc = sys->wrapDiff(pos[l] - pos[k]);
+ const Vector3 ac = bc+ab;
+ return compute_angle( ab.length2(), bc.length2(), ac.length2() );
+ }
+ HOST DEVICE inline float compute_angle(float distab2, float distbc2, float distac2) const {
// Find the cosine of the angle we want - <ABC
- float cos = (distab + distbc - distac2);
+ float cos = (distab2 + distbc2 - distac2);
- distab = 1.0f/sqrt(distab); //TODO: test other functions
- distbc = 1.0f/sqrt(distbc);
- cos *= 0.5f * distbc * distab;
+ distab2 = 1.0f/sqrt(distab2); //TODO: test other functions
+ distbc2 = 1.0f/sqrt(distbc2);
+ cos *= 0.5f * distbc2 * distab2;
// If the cosine is illegitimate, set it to 1 or -1 so that acos won't fail
if (cos < -1.0f) cos = -1.0f;
@@ -338,6 +347,9 @@ public:
else if (type == DIHEDRAL)
apply_dihedral_force(pos, sys, forces, particles[0], particles[1],
particles[2], particles[3], energy_deriv);
+ else if (type == VECANGLE)
+ apply_vecangle_force(pos, sys, forces, particles[0], particles[1],
+ particles[2], particles[3], energy_deriv);
}
__device__ inline void apply_bond_force(const Vector3* __restrict__ pos,
@@ -352,20 +364,18 @@ public:
#endif
}
- DEVICE inline void apply_angle_force(const Vector3* __restrict__ pos,
- const BaseGrid* __restrict__ sys,
- Vector3* __restrict__ force,
- int i, int j, int k, float energy_deriv) const {
-
-#ifdef __CUDA_ARCH__
- const Vector3 ab = -sys->wrapDiff(pos[j] - pos[i]);
- const Vector3 bc = -sys->wrapDiff(pos[k] - pos[j]);
- const Vector3 ac = sys->wrapDiff(pos[k] - pos[i]);
+ struct TwoVector3 {
+ Vector3 v1;
+ Vector3 v2;
+ };
+ DEVICE inline TwoVector3 get_angle_force(const Vector3& ab,
+ const Vector3& bc,
+ float energy_deriv) const {
// Find the distance between each pair of particles
float distab = ab.length2();
float distbc = bc.length2();
- const float distac2 = ac.length2();
+ const float distac2 = (ab+bc).length2();
// Find the cosine of the angle we want - <ABC
float cos = (distab + distbc - distac2);
@@ -382,12 +392,27 @@ public:
energy_deriv /= abs(sin) > 1e-3 ? sin : 1e-3; // avoid singularity
// Calculate the forces
- Vector3 force1 = -(energy_deriv*distab) * (ab * (cos*distab) + bc * distbc); // force on particle 1
- Vector3 force3 = (energy_deriv*distbc) * (bc * (cos*distbc) + ab * distab); // force on particle 3
+ TwoVector3 force;
+ force.v1 = (energy_deriv*distab) * (ab * (cos*distab) + bc * distbc); // force on 1st particle
+ force.v2 = -(energy_deriv*distbc) * (bc * (cos*distbc) + ab * distab); // force on last particle
+ return force;
+ }
- atomicAdd( &force[i], force1 );
- atomicAdd( &force[j], -(force1 + force3) );
- atomicAdd( &force[k], force3 );
+ DEVICE inline void apply_angle_force(const Vector3* __restrict__ pos,
+ const BaseGrid* __restrict__ sys,
+ Vector3* __restrict__ force,
+ int i, int j, int k, float energy_deriv) const {
+
+#ifdef __CUDA_ARCH__
+ const Vector3 ab = sys->wrapDiff(pos[j] - pos[i]);
+ const Vector3 bc = sys->wrapDiff(pos[k] - pos[j]);
+ // const Vector3 ac = sys->wrapDiff(pos[k] - pos[i]);
+
+ TwoVector3 f = get_angle_force(ab,bc, energy_deriv);
+
+ atomicAdd( &force[i], f.v1 );
+ atomicAdd( &force[j], -(f.v1 + f.v2) );
+ atomicAdd( &force[k], f.v2 );
#endif
}
@@ -431,6 +456,26 @@ public:
atomicAdd( &force[l], -f3 );
#endif
}
+ DEVICE inline void apply_vecangle_force(const Vector3* __restrict__ pos,
+ const BaseGrid* __restrict__ sys,
+ Vector3* __restrict__ force,
+ int i, int j, int k, int l, float energy_deriv) const {
+
+#ifdef __CUDA_ARCH__
+
+ const Vector3 ab = -sys->wrapDiff(pos[j] - pos[i]);
+ const Vector3 bc = -sys->wrapDiff(pos[k] - pos[j]);
+ // const Vector3 ac = sys->wrapDiff(pos[k] - pos[i]);
+
+ TwoVector3 f = get_angle_force(ab,bc, energy_deriv);
+
+ atomicAdd( &force[i], f.v1 );
+ atomicAdd( &force[j], -f.v1 );
+ atomicAdd( &force[k], -f.v2 );
+ atomicAdd( &force[l], f.v2 );
+#endif
+ }
+
};
#ifndef delgpuErrchk
--
GitLab