ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/UseTheForce/DarkSide/shapes.F90
(Generate patch)

Comparing trunk/OOPSE-4/src/UseTheForce/DarkSide/shapes.F90 (file contents):
Revision 1633 by gezelter, Fri Oct 22 20:22:48 2004 UTC vs.
Revision 1647 by chrisfen, Tue Oct 26 18:03:12 2004 UTC

# Line 25 | Line 25 | module shapes
25  
26    public :: do_shape_pair
27    public :: newShapeType
28 +  public :: complete_Shape_FF
29  
30  
31    type, private :: Shape
# Line 288 | Line 289 | contains  
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
# Line 378 | Line 326 | contains  
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        
# Line 388 | Line 334 | contains  
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)
# Line 1393 | Line 1339 | end subroutine makeShape
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 +

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines