| 35 |
|
* |
| 36 |
|
* [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). |
| 37 |
|
* [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). |
| 38 |
< |
* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). |
| 39 |
< |
* [4] Vardeman & Gezelter, in progress (2009). |
| 38 |
> |
* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008). |
| 39 |
> |
* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). |
| 40 |
> |
* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). |
| 41 |
|
*/ |
| 42 |
|
|
| 43 |
|
#include <stdio.h> |
| 86 |
|
// add it to the map: |
| 87 |
|
AtomTypeProperties atp = atomType->getATP(); |
| 88 |
|
|
| 89 |
< |
pair<map<int,ShapeAtomType*>::iterator, bool> ret; |
| 90 |
< |
ret = shapesMap.insert( pair<int, ShapeAtomType*>(atp.ident, atomType)); |
| 91 |
< |
if (ret.second == false) { |
| 92 |
< |
sprintf( painCave.errmsg, |
| 93 |
< |
"SHAPES already had a previous entry with ident %d\n", |
| 94 |
< |
atp.ident); |
| 95 |
< |
painCave.severity = OPENMD_INFO; |
| 96 |
< |
painCave.isFatal = 0; |
| 97 |
< |
simError(); |
| 98 |
< |
} |
| 99 |
< |
|
| 100 |
< |
ShapesMap.insert( pair<int, ShapeAtomType*>(atp.ident, sAtomType) ); |
| 101 |
< |
|
| 102 |
< |
} else if (atomType->isLennardJones()) { |
| 103 |
< |
d1 = getLJSigma(atomType) / sqrt(2.0); |
| 104 |
< |
e1 = getLJEpsilon(atomType); |
| 89 |
> |
if (atomType->isShape() ) { |
| 90 |
> |
pair<map<int,ShapeAtomType*>::iterator, bool> ret; |
| 91 |
> |
ret = ShapesMap.insert( pair<int, ShapeAtomType*>(atp.ident, atomType)); |
| 92 |
> |
if (ret.second == false) { |
| 93 |
> |
sprintf( painCave.errMsg, |
| 94 |
> |
"SHAPES already had a previous entry with ident %d\n", |
| 95 |
> |
atp.ident); |
| 96 |
> |
painCave.severity = OPENMD_INFO; |
| 97 |
> |
painCave.isFatal = 0; |
| 98 |
> |
simError(); |
| 99 |
> |
} |
| 100 |
> |
|
| 101 |
> |
ShapesMap.insert( pair<int, ShapeAtomType*>(atp.ident, static_cast<ShapeAtomType*>(atomType)) ); |
| 102 |
> |
|
| 103 |
> |
} else if (atomType->isLennardJones()) { |
| 104 |
> |
RealType d1 = getLJSigma(atomType) / sqrt(2.0); |
| 105 |
> |
RealType e1 = getLJEpsilon(atomType); |
| 106 |
|
} else { |
| 107 |
|
sprintf( painCave.errMsg, |
| 108 |
|
"SHAPES::addType was passed an atomType (%s) that does not\n" |
| 263 |
|
// sin(theta) is near 0.0: |
| 264 |
|
|
| 265 |
|
RealType sti2 = 1.0 - cti*cti; |
| 266 |
+ |
RealType proji; |
| 267 |
+ |
Vector3d dcpidr, dcpidu, dspidr, dspidu; |
| 268 |
|
if (fabs(sti2) < 1.0e-12) { |
| 269 |
< |
RealType proji = sqrt(r * 1.0e-12); |
| 270 |
< |
Vector3d dcpidx(1.0 / proji, |
| 271 |
< |
0.0, |
| 272 |
< |
|
| 273 |
< |
dcpidx = 1.0_dp / proji |
| 274 |
< |
dcpidy = 0.0_dp |
| 275 |
< |
dcpidux = xi / proji |
| 276 |
< |
dcpiduy = 0.0_dp |
| 277 |
< |
dspidx = 0.0_dp |
| 278 |
< |
dspidy = 1.0_dp / proji |
| 279 |
< |
dspidux = 0.0_dp |
| 280 |
< |
dspiduy = yi / proji |
| 281 |
< |
else |
| 282 |
< |
proji = sqrt(xi2 + yi2) |
| 283 |
< |
proji3 = proji*proji*proji |
| 284 |
< |
dcpidx = 1.0_dp / proji - xi2 / proji3 |
| 285 |
< |
dcpidy = - xi * yi / proji3 |
| 286 |
< |
dcpidux = xi / proji - (xi2 * xi) / proji3 |
| 287 |
< |
dcpiduy = - (xi * yi2) / proji3 |
| 288 |
< |
dspidx = - xi * yi / proji3 |
| 289 |
< |
dspidy = 1.0_dp / proji - yi2 / proji3 |
| 286 |
< |
dspidux = - (yi * xi2) / proji3 |
| 287 |
< |
dspiduy = yi / proji - (yi2 * yi) / proji3 |
| 288 |
< |
endif |
| 269 |
> |
proji = sqrt(r * 1.0e-12); |
| 270 |
> |
dcpidr = Vector3d( 1.0 / proji, 0.0, 0.0); |
| 271 |
> |
dcpidu = Vector3d( xi / proji, 0.0, 0.0); |
| 272 |
> |
dspidr = Vector3d( 0.0, 1.0 / proji, 0.0); |
| 273 |
> |
dspidu = Vector3d (0.0, yi / proji, 0.0); |
| 274 |
> |
} else { |
| 275 |
> |
proji = sqrt(xi2 + yi2); |
| 276 |
> |
RealType proji3 = proji*proji*proji; |
| 277 |
> |
dcpidr = Vector3d(1.0_dp / proji - xi2 / proji3, |
| 278 |
> |
- xi * yi / proji3, |
| 279 |
> |
0.0); |
| 280 |
> |
dcpidu = Vector3d( xi / proji - (xi2 * xi) / proji3, |
| 281 |
> |
- (xi * yi2) / proji3, |
| 282 |
> |
0.0); |
| 283 |
> |
dspidr = Vector3d(- xi * yi / proji3, |
| 284 |
> |
1.0_dp / proji - yi2 / proji3, |
| 285 |
> |
0.0); |
| 286 |
> |
dspidu = Vector3d(- (yi * xi2) / proji3, |
| 287 |
> |
yi / proji - (yi2 * yi) / proji3, |
| 288 |
> |
0.0); |
| 289 |
> |
} |
| 290 |
|
|
| 291 |
< |
cpi = xi / proji |
| 292 |
< |
dcpidz = 0.0_dp |
| 293 |
< |
dcpiduz = 0.0_dp |
| 291 |
> |
cpi = xi / proji; |
| 292 |
> |
dcpidr.z() = 0.0; |
| 293 |
> |
dcpidu.z() = 0.0; |
| 294 |
> |
|
| 295 |
> |
spi = yi / proji; |
| 296 |
> |
dspidr.z() = 0.0; |
| 297 |
> |
dspidu.z() = 0.0; |
| 298 |
|
|
| 294 |
– |
spi = yi / proji |
| 295 |
– |
dspidz = 0.0_dp |
| 296 |
– |
dspiduz = 0.0_dp |
| 299 |
|
|
| 298 |
– |
|
| 299 |
– |
|
| 300 |
– |
|
| 300 |
|
RealType sigma0 = mixer.sigma0; |
| 301 |
|
RealType dw = mixer.dw; |
| 302 |
|
RealType eps0 = mixer.eps0; |
