25 |
|
|
26 |
|
public :: do_shape_pair |
27 |
|
public :: newShapeType |
28 |
+ |
public :: complete_Shape_FF |
29 |
|
|
30 |
|
|
31 |
|
type, private :: Shape |
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 |
295 |
< |
real(kind=dp) :: costheta, cpi, theta, Pi, phi, thisDP, sigma |
295 |
> |
real(kind=dp) :: thisDP, sigma |
296 |
|
integer :: alloc_stat, iTheta, iPhi, nSteps, nAtypes, thisIP, current |
297 |
|
logical :: thisProperty |
298 |
|
|
299 |
– |
Pi = 4.0d0 * datan(1.0d0) |
300 |
– |
|
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 |
< |
|
308 |
< |
bigSigma = 0.0d0 |
309 |
< |
do i = 1, ShapeMap%currentShape |
310 |
< |
|
311 |
< |
! Scan over theta and phi to |
312 |
< |
! find the largest contact in any direction.... |
313 |
< |
|
314 |
< |
myBigSigma = 0.0d0 |
315 |
< |
|
316 |
< |
do iTheta = 0, nSteps |
317 |
< |
theta = (Pi/2.0d0)*(dble(iTheta)/dble(nSteps)) |
318 |
< |
costheta = cos(theta) |
319 |
< |
|
320 |
< |
call Associated_Legendre(costheta, ShapeMap%Shapes(i)%bigL, & |
321 |
< |
ShapeMap%Shapes(i)%bigM, lmax, plm_i, dlm_i) |
322 |
< |
|
323 |
< |
do iPhi = 0, nSteps |
324 |
< |
phi = -Pi + 2.0d0 * Pi * (dble(iPhi)/dble(nSteps)) |
325 |
< |
cpi = cos(phi) |
326 |
< |
spi = sin(phi) |
327 |
< |
|
328 |
< |
call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(i)%bigM, & |
329 |
< |
CHEBYSHEV_TN, tm_i, dtm_i) |
330 |
< |
call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(i)%bigM, & |
331 |
< |
CHEBYSHEV_UN, um_i, dum_i) |
332 |
< |
|
333 |
< |
thisSigma = 0.0d0 |
334 |
< |
|
335 |
< |
do lm = 1, ShapeMap%Shapes(i)%nContactFuncs |
336 |
< |
|
337 |
< |
l = ShapeMap%Shapes(i)%ContactFuncLValue(lm) |
338 |
< |
m = ShapeMap%Shapes(i)%ContactFuncMValue(lm) |
339 |
< |
coeff = ShapeMap%Shapes(i)%ContactFuncCoefficient(lm) |
340 |
< |
function_type = ShapeMap%Shapes(i)%ContactFunctionType(lm) |
341 |
< |
|
342 |
< |
if ((function_type .eq. SH_COS).or.(m.eq.0)) then |
343 |
< |
Phunc = coeff * tm_i(m) |
344 |
< |
else |
345 |
< |
Phunc = coeff * spi * um_i(m-1) |
346 |
< |
endif |
347 |
< |
|
348 |
< |
thisSigma = thisSigma + plm_i(l,m)*Phunc |
349 |
< |
enddo |
350 |
< |
|
351 |
< |
if (thisSigma.gt.myBigSigma) myBigSigma = thisSigma |
352 |
< |
enddo |
353 |
< |
enddo |
354 |
< |
|
355 |
< |
if (myBigSigma.gt.bigSigma) bigSigma = myBigSigma |
356 |
< |
enddo |
357 |
< |
|
305 |
> |
|
306 |
|
nAtypes = getSize(atypes) |
307 |
|
|
308 |
|
if (nAtypes == 0) then |
326 |
|
ShapeMap%Shapes(current)%isLJ = .true. |
327 |
|
|
328 |
|
ShapeMap%Shapes(current)%epsilon = getEpsilon(thisIP) |
329 |
< |
sigma = getSigma(thisIP) |
382 |
< |
ShapeMap%Shapes(current)%sigma = sigma |
383 |
< |
if (sigma .gt. bigSigma) bigSigma = thisDP |
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 |
+ |
|