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 |
+ |
|