| 6 |
|
use vector_class |
| 7 |
|
use simulation |
| 8 |
|
use status |
| 9 |
+ |
use lj |
| 10 |
|
#ifdef IS_MPI |
| 11 |
|
use mpiSimulation |
| 12 |
|
#endif |
| 25 |
|
|
| 26 |
|
public :: do_shape_pair |
| 27 |
|
public :: newShapeType |
| 28 |
+ |
public :: complete_Shape_FF |
| 29 |
|
|
| 30 |
|
|
| 31 |
|
type, private :: Shape |
| 106 |
|
call getMatchingElementList(atypes, "is_Shape", .true., & |
| 107 |
|
nShapeTypes, MatchList) |
| 108 |
|
|
| 109 |
< |
call getMatchingElementList(atypes, "is_LJ", .true., & |
| 109 |
> |
call getMatchingElementList(atypes, "is_LennardJones", .true., & |
| 110 |
|
nLJTypes, MatchList) |
| 111 |
|
|
| 112 |
|
ShapeMap%n_shapes = nShapeTypes + nLJTypes |
| 289 |
|
|
| 290 |
|
end subroutine allocateShape |
| 291 |
|
|
| 292 |
< |
subroutine init_Shape_FF(status) |
| 292 |
> |
subroutine complete_Shape_FF(status) |
| 293 |
|
integer :: status |
| 294 |
|
integer :: i, j, l, m, lm, function_type |
| 295 |
< |
real(kind=dp) :: bigSigma, myBigSigma, thisSigma, coeff, Phunc, spi |
| 294 |
< |
real(kind=dp) :: costheta, cpi, theta, Pi, phi, thisDP |
| 295 |
> |
real(kind=dp) :: thisDP, sigma |
| 296 |
|
integer :: alloc_stat, iTheta, iPhi, nSteps, nAtypes, thisIP, current |
| 297 |
|
logical :: thisProperty |
| 298 |
|
|
| 298 |
– |
Pi = 4.0d0 * datan(1.0d0) |
| 299 |
– |
|
| 299 |
|
status = 0 |
| 300 |
|
if (ShapeMap%currentShape == 0) then |
| 301 |
|
call handleError("init_Shape_FF", "No members in ShapeMap") |
| 302 |
|
status = -1 |
| 303 |
|
return |
| 304 |
|
end if |
| 305 |
< |
|
| 307 |
< |
bigSigma = 0.0d0 |
| 308 |
< |
do i = 1, ShapeMap%currentShape |
| 309 |
< |
|
| 310 |
< |
! Scan over theta and phi to |
| 311 |
< |
! find the largest contact in any direction.... |
| 312 |
< |
|
| 313 |
< |
myBigSigma = 0.0d0 |
| 314 |
< |
|
| 315 |
< |
do iTheta = 0, nSteps |
| 316 |
< |
theta = (Pi/2.0d0)*(dble(iTheta)/dble(nSteps)) |
| 317 |
< |
costheta = cos(theta) |
| 318 |
< |
|
| 319 |
< |
call Associated_Legendre(costheta, ShapeMap%Shapes(i)%bigL, & |
| 320 |
< |
ShapeMap%Shapes(i)%bigM, lmax, plm_i, dlm_i) |
| 321 |
< |
|
| 322 |
< |
do iPhi = 0, nSteps |
| 323 |
< |
phi = -Pi + 2.0d0 * Pi * (dble(iPhi)/dble(nSteps)) |
| 324 |
< |
cpi = cos(phi) |
| 325 |
< |
spi = sin(phi) |
| 326 |
< |
|
| 327 |
< |
call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(i)%bigM, & |
| 328 |
< |
CHEBYSHEV_TN, tm_i, dtm_i) |
| 329 |
< |
call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(i)%bigM, & |
| 330 |
< |
CHEBYSHEV_UN, um_i, dum_i) |
| 331 |
< |
|
| 332 |
< |
thisSigma = 0.0d0 |
| 333 |
< |
|
| 334 |
< |
do lm = 1, ShapeMap%Shapes(i)%nContactFuncs |
| 335 |
< |
|
| 336 |
< |
l = ShapeMap%Shapes(i)%ContactFuncLValue(lm) |
| 337 |
< |
m = ShapeMap%Shapes(i)%ContactFuncMValue(lm) |
| 338 |
< |
coeff = ShapeMap%Shapes(i)%ContactFuncCoefficient(lm) |
| 339 |
< |
function_type = ShapeMap%Shapes(i)%ContactFunctionType(lm) |
| 340 |
< |
|
| 341 |
< |
if ((function_type .eq. SH_COS).or.(m.eq.0)) then |
| 342 |
< |
Phunc = coeff * tm_i(m) |
| 343 |
< |
else |
| 344 |
< |
Phunc = coeff * spi * um_i(m-1) |
| 345 |
< |
endif |
| 346 |
< |
|
| 347 |
< |
thisSigma = thisSigma + plm_i(l,m)*Phunc |
| 348 |
< |
enddo |
| 349 |
< |
|
| 350 |
< |
if (thisSigma.gt.myBigSigma) myBigSigma = thisSigma |
| 351 |
< |
enddo |
| 352 |
< |
enddo |
| 353 |
< |
|
| 354 |
< |
if (myBigSigma.gt.bigSigma) bigSigma = myBigSigma |
| 355 |
< |
enddo |
| 356 |
< |
|
| 305 |
> |
|
| 306 |
|
nAtypes = getSize(atypes) |
| 307 |
|
|
| 308 |
|
if (nAtypes == 0) then |
| 312 |
|
|
| 313 |
|
do i = 1, nAtypes |
| 314 |
|
|
| 315 |
< |
call getElementProperty(atypes, i, "is_LJ", thisProperty) |
| 315 |
> |
call getElementProperty(atypes, i, "is_LennardJones", thisProperty) |
| 316 |
|
|
| 317 |
|
if (thisProperty) then |
| 318 |
|
|
| 325 |
|
|
| 326 |
|
ShapeMap%Shapes(current)%isLJ = .true. |
| 327 |
|
|
| 328 |
< |
call getElementProperty(atypes, i, "lj_epsilon", thisDP) |
| 329 |
< |
ShapeMap%Shapes(current)%epsilon = thisDP |
| 381 |
< |
|
| 382 |
< |
call getElementProperty(atypes, i, "lj_sigma", thisDP) |
| 383 |
< |
ShapeMap%Shapes(current)%sigma = thisDP |
| 384 |
< |
if (thisDP .gt. bigSigma) bigSigma = thisDP |
| 328 |
> |
ShapeMap%Shapes(current)%epsilon = getEpsilon(thisIP) |
| 329 |
> |
ShapeMap%Shapes(current)%sigma = getSigma(thisIP) |
| 330 |
|
|
| 331 |
|
endif |
| 332 |
|
|
| 334 |
|
|
| 335 |
|
haveShapeMap = .true. |
| 336 |
|
|
| 337 |
< |
end subroutine init_Shape_FF |
| 337 |
> |
end subroutine complete_Shape_FF |
| 338 |
|
|
| 339 |
|
subroutine do_shape_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, & |
| 340 |
|
pot, A, f, t, do_pot) |
| 1339 |
|
|
| 1340 |
|
return |
| 1341 |
|
end subroutine makeShape |
| 1342 |
+ |
|
| 1343 |
+ |
subroutine completeShapeFF(status) |
| 1344 |
+ |
|
| 1345 |
+ |
use shapes, only: complete_Shape_FF |
| 1346 |
+ |
|
| 1347 |
+ |
integer, intent(out) :: status |
| 1348 |
+ |
integer :: myStatus |
| 1349 |
+ |
|
| 1350 |
+ |
myStatus = 0 |
| 1351 |
+ |
|
| 1352 |
+ |
call complete_Shape_FF(myStatus) |
| 1353 |
+ |
|
| 1354 |
+ |
status = myStatus |
| 1355 |
+ |
|
| 1356 |
+ |
return |
| 1357 |
+ |
end subroutine completeShapeFF |
| 1358 |
+ |
|
