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 1608 by gezelter, Wed Oct 20 04:02:48 2004 UTC vs.
Revision 2211 by chrisfen, Thu Apr 21 14:12:19 2005 UTC

# Line 1 | Line 1
1 + !!
2 + !! Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3 + !!
4 + !! The University of Notre Dame grants you ("Licensee") a
5 + !! non-exclusive, royalty free, license to use, modify and
6 + !! redistribute this software in source and binary code form, provided
7 + !! that the following conditions are met:
8 + !!
9 + !! 1. Acknowledgement of the program authors must be made in any
10 + !!    publication of scientific results based in part on use of the
11 + !!    program.  An acceptable form of acknowledgement is citation of
12 + !!    the article in which the program was described (Matthew
13 + !!    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 + !!    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 + !!    Parallel Simulation Engine for Molecular Dynamics,"
16 + !!    J. Comput. Chem. 26, pp. 252-271 (2005))
17 + !!
18 + !! 2. Redistributions of source code must retain the above copyright
19 + !!    notice, this list of conditions and the following disclaimer.
20 + !!
21 + !! 3. Redistributions in binary form must reproduce the above copyright
22 + !!    notice, this list of conditions and the following disclaimer in the
23 + !!    documentation and/or other materials provided with the
24 + !!    distribution.
25 + !!
26 + !! This software is provided "AS IS," without a warranty of any
27 + !! kind. All express or implied conditions, representations and
28 + !! warranties, including any implied warranty of merchantability,
29 + !! fitness for a particular purpose or non-infringement, are hereby
30 + !! excluded.  The University of Notre Dame and its licensors shall not
31 + !! be liable for any damages suffered by licensee as a result of
32 + !! using, modifying or distributing the software or its
33 + !! derivatives. In no event will the University of Notre Dame or its
34 + !! licensors be liable for any lost revenue, profit or data, or for
35 + !! direct, indirect, special, consequential, incidental or punitive
36 + !! damages, however caused and regardless of the theory of liability,
37 + !! arising out of the use of or inability to use software, even if the
38 + !! University of Notre Dame has been advised of the possibility of
39 + !! such damages.
40 + !!
41 +
42 +
43   module shapes
44  
45    use force_globals
# Line 6 | Line 48 | module shapes
48    use vector_class
49    use simulation
50    use status
51 +  use lj
52   #ifdef IS_MPI
53    use mpiSimulation
54   #endif
55    implicit none
56  
57    PRIVATE
58 <  
58 >
59    INTEGER, PARAMETER:: CHEBYSHEV_TN = 1
60    INTEGER, PARAMETER:: CHEBYSHEV_UN = 2
61    INTEGER, PARAMETER:: LAGUERRE     = 3
# Line 24 | Line 67 | module shapes
67  
68    public :: do_shape_pair
69    public :: newShapeType
70 +  public :: complete_Shape_FF
71 +  public :: destroyShapeTypes
72  
28 –
73    type, private :: Shape
74       integer :: atid
75       integer :: nContactFuncs
# Line 49 | Line 93 | module shapes
93       real ( kind = dp )  :: epsilon
94       real ( kind = dp )  :: sigma
95    end type Shape
96 <  
96 >
97    type, private :: ShapeList
98       integer :: n_shapes = 0
99       integer :: currentShape = 0
100       type (Shape), pointer :: Shapes(:)      => null()
101       integer, pointer      :: atidToShape(:) => null()
102    end type ShapeList
103 <  
103 >
104    type(ShapeList), save :: ShapeMap
105  
106    integer :: lmax
63 –  real (kind=dp), allocatable, dimension(:,:) :: plm_i, dlm_i, plm_j, dlm_j
64 –  real (kind=dp), allocatable, dimension(:) :: tm_i, dtm_i, um_i, dum_i
65 –  real (kind=dp), allocatable, dimension(:) :: tm_j, dtm_j, um_j, dum_j
107  
108   contains  
109  
# Line 71 | Line 112 | contains  
112         nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
113         RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
114         StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
115 <       myAtid, status)
115 >       c_ident, status)
116  
117      integer :: nContactFuncs
118      integer :: nRangeFuncs
119      integer :: nStrengthFuncs
120      integer :: shape_ident
121      integer :: status
122 <    integer :: myAtid
122 >    integer :: c_ident
123 >    integer :: myATID
124      integer :: bigL
125      integer :: bigM
126      integer :: j, me, nShapeTypes, nLJTypes, ntypes, current, alloc_stat
# Line 103 | Line 145 | contains  
145  
146         call getMatchingElementList(atypes, "is_Shape", .true., &
147              nShapeTypes, MatchList)
148 <      
149 <       call getMatchingElementList(atypes, "is_LJ", .true., &
148 >
149 >       call getMatchingElementList(atypes, "is_LennardJones", .true., &
150              nLJTypes, MatchList)
151 <      
151 >
152         ShapeMap%n_shapes = nShapeTypes + nLJTypes
153 <      
153 >
154         allocate(ShapeMap%Shapes(nShapeTypes + nLJTypes))
155 <      
155 >
156         ntypes = getSize(atypes)
157 <      
157 >
158         allocate(ShapeMap%atidToShape(ntypes))
159      end if
160 <    
160 >
161      ShapeMap%currentShape = ShapeMap%currentShape + 1
162      current = ShapeMap%currentShape
163  
# Line 125 | Line 167 | contains  
167         status = -1
168         return
169      endif
170 +    
171 +    myATID = getFirstMatchingElement(atypes, 'c_ident', c_ident)
172 + !    write(*,*) 'myATID= ', myATID, ' c_ident = ', c_ident
173  
174 <    call getElementProperty(atypes, myAtid, "c_ident", me)
175 <    ShapeMap%atidToShape(me)                         = current
131 <    ShapeMap%Shapes(current)%atid                    = me
174 >    ShapeMap%atidToShape(myATID)                     = current
175 >    ShapeMap%Shapes(current)%atid                    = myATID
176      ShapeMap%Shapes(current)%nContactFuncs           = nContactFuncs
177      ShapeMap%Shapes(current)%nRangeFuncs             = nRangeFuncs
178      ShapeMap%Shapes(current)%nStrengthFuncs          = nStrengthFuncs
# Line 147 | Line 191 | contains  
191  
192      bigL = -1
193      bigM = -1
194 <    
194 >
195      do j = 1, ShapeMap%Shapes(current)%nContactFuncs
196         if (ShapeMap%Shapes(current)%ContactFuncLValue(j) .gt. bigL) then
197            bigL = ShapeMap%Shapes(current)%ContactFuncLValue(j)
# Line 185 | Line 229 | contains  
229      type(Shape), intent(inout) :: myShape
230      integer, intent(out) :: stat
231      integer :: alloc_stat
232 <
232 >
233 >    stat = 0
234      if (associated(myShape%contactFuncLValue)) then
235         deallocate(myShape%contactFuncLValue)
236      endif
# Line 251 | Line 296 | contains  
296         stat = -1
297         return
298      endif
299 <    
299 >
300      if (associated(myShape%strengthFuncLValue)) then
301         deallocate(myShape%strengthFuncLValue)
302      endif
# Line 285 | Line 330 | contains  
330         return
331      endif
332  
333 +    return
334 +
335    end subroutine allocateShape
336 <    
337 <  subroutine init_Shape_FF(status)
336 >
337 >  subroutine complete_Shape_FF(status)
338      integer :: status
339      integer :: i, j, l, m, lm, function_type
340 <    real(kind=dp) :: bigSigma, myBigSigma, thisSigma, coeff, Phunc, spi
341 <    real(kind=dp) :: costheta, cpi, theta, Pi, phi, thisDP
295 <    integer :: alloc_stat, iTheta, iPhi, nSteps, nAtypes, thisIP, current
340 >    real(kind=dp) :: thisDP, sigma
341 >    integer :: alloc_stat, iTheta, iPhi, nSteps, nAtypes, myATID, current
342      logical :: thisProperty
343  
298 –    Pi = 4.0d0 * datan(1.0d0)
299 –
344      status = 0
345      if (ShapeMap%currentShape == 0) then
346         call handleError("init_Shape_FF", "No members in ShapeMap")
347         status = -1
348         return
349      end if
306 –
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)
350  
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 –
351      nAtypes = getSize(atypes)
352  
353      if (nAtypes == 0) then
354         status = -1
355         return
356 <    end if
356 >    end if      
357  
358 +    ! atypes comes from c side
359      do i = 1, nAtypes
360 +    
361 +       myATID = getFirstMatchingElement(atypes, 'c_ident', i)
362 +       call getElementProperty(atypes, myATID, "is_LennardJones", thisProperty)
363 + !       write(*,*) 'is_LJ = ', thisProperty, ' for atid = ', myATID
364        
366 –       call getElementProperty(atypes, i, "is_LJ", thisProperty)
367 –      
365         if (thisProperty) then
366 <          
366 >                
367            ShapeMap%currentShape = ShapeMap%currentShape + 1
368            current = ShapeMap%currentShape
369  
370 <          call getElementProperty(atypes, i, "c_ident",  thisIP)
371 <          ShapeMap%atidToShape(thisIP) = current
375 <          ShapeMap%Shapes(current)%atid = thisIP
370 >          ShapeMap%atidToShape(myATID) = current
371 >          ShapeMap%Shapes(current)%atid = myATID
372  
373            ShapeMap%Shapes(current)%isLJ = .true.
374  
375 <          call getElementProperty(atypes, i, "lj_epsilon", thisDP)
376 <          ShapeMap%Shapes(current)%epsilon = thisDP
375 >          ShapeMap%Shapes(current)%epsilon = getEpsilon(myATID)
376 >          ShapeMap%Shapes(current)%sigma = getSigma(myATID)
377  
382 –          call getElementProperty(atypes, i, "lj_sigma",   thisDP)
383 –          ShapeMap%Shapes(current)%sigma = thisDP          
384 –          if (thisDP .gt. bigSigma) bigSigma = thisDP
385 –          
378         endif
379 <      
379 >
380      end do
381  
382      haveShapeMap = .true.
383 <    
384 <  end subroutine init_Shape_FF
385 <    
383 >
384 > !    do i = 1, ShapeMap%n_shapes
385 > !       write(*,*) 'i = ', i, ' isLJ = ', ShapeMap%Shapes(i)%isLJ
386 > !    end do
387 >
388 >  end subroutine complete_Shape_FF
389 >
390    subroutine do_shape_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
391         pot, A, f, t, do_pot)
392 <    
392 >
393 >    INTEGER, PARAMETER:: LMAX         = 64
394 >    INTEGER, PARAMETER:: MMAX         = 64
395 >
396      integer, intent(in) :: atom1, atom2
397      real (kind=dp), intent(inout) :: rij, r2
398      real (kind=dp), dimension(3), intent(in) :: d
399      real (kind=dp), dimension(3), intent(inout) :: fpair
400 <    real (kind=dp) :: pot, vpair, sw
400 >    real (kind=dp) :: pot, vpair, sw, dswdr
401      real (kind=dp), dimension(9,nLocal) :: A
402      real (kind=dp), dimension(3,nLocal) :: f
403      real (kind=dp), dimension(3,nLocal) :: t
# Line 409 | Line 408 | contains  
408      integer :: l, m, lm, id1, id2, localError, function_type
409      real (kind=dp) :: sigma_i, s_i, eps_i, sigma_j, s_j, eps_j
410      real (kind=dp) :: coeff
411 +    real (kind=dp) :: pot_temp
412  
413      real (kind=dp) :: dsigmaidx, dsigmaidy, dsigmaidz
414      real (kind=dp) :: dsigmaidux, dsigmaiduy, dsigmaiduz
# Line 427 | Line 427 | contains  
427  
428      real (kind=dp) :: xi, yi, zi, xj, yj, zj, xi2, yi2, zi2, xj2, yj2, zj2
429  
430 +    real (kind=dp) :: sti2, stj2
431 +
432      real (kind=dp) :: proji, proji3, projj, projj3
433      real (kind=dp) :: cti, ctj, cpi, cpj, spi, spj
434      real (kind=dp) :: Phunc, sigma, s, eps, rtdenom, rt
# Line 458 | Line 460 | contains  
460      real (kind=dp) :: dsduxi, dsduyi, dsduzi
461      real (kind=dp) :: dsdxj, dsdyj, dsdzj
462      real (kind=dp) :: dsduxj, dsduyj, dsduzj
463 <    
463 >
464      real (kind=dp) :: depsdxi, depsdyi, depsdzi
465      real (kind=dp) :: depsduxi, depsduyi, depsduzi
466      real (kind=dp) :: depsdxj, depsdyj, depsdzj
# Line 485 | Line 487 | contains  
487      real (kind=dp) :: fxji, fyji, fzji, fxjj, fyjj, fzjj
488      real (kind=dp) :: fxradial, fyradial, fzradial
489  
490 +    real (kind=dp) :: plm_i(0:LMAX,0:MMAX), dlm_i(0:LMAX,0:MMAX)
491 +    real (kind=dp) :: plm_j(0:LMAX,0:MMAX), dlm_j(0:LMAX,0:MMAX)
492 +    real (kind=dp) :: tm_i(0:MMAX), dtm_i(0:MMAX), um_i(0:MMAX), dum_i(0:MMAX)
493 +    real (kind=dp) :: tm_j(0:MMAX), dtm_j(0:MMAX), um_j(0:MMAX), dum_j(0:MMAX)
494 +
495      if (.not.haveShapeMap) then
496         call handleError("calc_shape", "NO SHAPEMAP!!!!")
497         return      
498      endif
499 <    
499 >
500      !! We assume that the rotation matrices have already been calculated
501      !! and placed in the A array.
495 –        
502      r3 = r2*rij
503      r5 = r3*r2
504 <    
504 >
505      drdxi = -d(1) / rij
506      drdyi = -d(2) / rij
507      drdzi = -d(3) / rij
# Line 503 | Line 509 | contains  
509      drdxj = d(1) / rij
510      drdyj = d(2) / rij
511      drdzj = d(3) / rij
512 <    
512 >
513      ! find the atom type id (atid) for each atom:
514   #ifdef IS_MPI
515      atid1 = atid_Row(atom1)
# Line 514 | Line 520 | contains  
520   #endif
521  
522      ! use the atid to find the shape type (st) for each atom:
517 –
523      st1 = ShapeMap%atidToShape(atid1)
524      st2 = ShapeMap%atidToShape(atid2)
525      
526 + !    write(*,*) atom1, atom2, atid1, atid2, st1, st2, ShapeMap%Shapes(st1)%isLJ, ShapeMap%Shapes(st2)%isLJ
527 +
528      if (ShapeMap%Shapes(st1)%isLJ) then
529 +
530         sigma_i = ShapeMap%Shapes(st1)%sigma
531         s_i = ShapeMap%Shapes(st1)%sigma
532         eps_i = ShapeMap%Shapes(st1)%epsilon
# Line 545 | Line 553 | contains  
553   #ifdef IS_MPI
554         ! rotate the inter-particle separation into the two different
555         ! body-fixed coordinate systems:
556 <      
556 >
557         xi = A_row(1,atom1)*d(1) + A_row(2,atom1)*d(2) + A_row(3,atom1)*d(3)
558         yi = A_row(4,atom1)*d(1) + A_row(5,atom1)*d(2) + A_row(6,atom1)*d(3)
559         zi = A_row(7,atom1)*d(1) + A_row(8,atom1)*d(2) + A_row(9,atom1)*d(3)
560 <      
560 >
561   #else
562         ! rotate the inter-particle separation into the two different
563         ! body-fixed coordinate systems:
564 <      
564 >
565         xi = a(1,atom1)*d(1) + a(2,atom1)*d(2) + a(3,atom1)*d(3)
566         yi = a(4,atom1)*d(1) + a(5,atom1)*d(2) + a(6,atom1)*d(3)
567         zi = a(7,atom1)*d(1) + a(8,atom1)*d(2) + a(9,atom1)*d(3)
568 <      
568 >
569   #endif
570 <      
570 >
571         xi2 = xi*xi
572         yi2 = yi*yi
573 <       zi2 = zi*zi
566 <      
567 <       proji = sqrt(xi2 + yi2)
568 <       proji3 = proji*proji*proji
569 <      
573 >       zi2 = zi*zi            
574         cti = zi / rij
575 +
576 +       if (cti .gt. 1.0_dp) cti = 1.0_dp
577 +       if (cti .lt. -1.0_dp) cti = -1.0_dp
578 +
579         dctidx = - zi * xi / r3
580         dctidy = - zi * yi / r3
581         dctidz = 1.0d0 / rij - zi2 / r3
582 <       dctidux =  yi / rij
583 <       dctiduy = -xi / rij
584 <       dctiduz = 0.0d0
585 <      
582 >       dctidux = - (zi * xi2) / r3
583 >       dctiduy = - (zi * yi2) / r3
584 >       dctiduz = zi / rij - (zi2 * zi) / r3
585 >
586 >       ! this is an attempt to try to truncate the singularity when
587 >       ! sin(theta) is near 0.0:
588 >
589 >       sti2 = 1.0_dp - cti*cti
590 >       if (dabs(sti2) .lt. 1.0d-12) then
591 >          proji = sqrt(rij * 1.0d-12)
592 >          dcpidx = 1.0d0 / proji
593 >          dcpidy = 0.0d0
594 >          dcpidux = xi / proji
595 >          dcpiduy = 0.0d0
596 >          dspidx = 0.0d0
597 >          dspidy = 1.0d0 / proji
598 >          dspidux = 0.0d0
599 >          dspiduy = yi / proji
600 >       else
601 >          proji = sqrt(xi2 + yi2)
602 >          proji3 = proji*proji*proji
603 >          dcpidx = 1.0d0 / proji - xi2 / proji3
604 >          dcpidy = - xi * yi / proji3
605 >          dcpidux = xi / proji - (xi2 * xi) / proji3
606 >          dcpiduy = - (xi * yi2) / proji3
607 >          dspidx = - xi * yi / proji3
608 >          dspidy = 1.0d0 / proji - yi2 / proji3
609 >          dspidux = - (yi * xi2) / proji3
610 >          dspiduy = yi / proji - (yi2 * yi) / proji3
611 >       endif
612 >
613         cpi = xi / proji
579 –       dcpidx = 1.0d0 / proji - xi2 / proji3
580 –       dcpidy = - xi * yi / proji3
614         dcpidz = 0.0d0
615 <       dcpidux = xi * yi * zi / proji3
616 <       dcpiduy = -zi * (1.0d0 / proji - xi2 / proji3)
584 <       dcpiduz = -yi * (1.0d0 / proji - xi2 / proji3)  - (xi2 * yi / proji3)
585 <      
615 >       dcpiduz = 0.0d0
616 >
617         spi = yi / proji
587 –       dspidx = - xi * yi / proji3
588 –       dspidy = 1.0d0 / proji - yi2 / proji3
618         dspidz = 0.0d0
619 <       dspidux = -zi * (1.0d0 / proji - yi2 / proji3)
591 <       dspiduy = xi * yi * zi / proji3
592 <       dspiduz = xi * (1.0d0 / proji - yi2 / proji3) + (xi * yi2 / proji3)
619 >       dspiduz = 0.0d0
620  
621 <       call Associated_Legendre(cti, ShapeMap%Shapes(st1)%bigL, &
622 <            ShapeMap%Shapes(st1)%bigM, lmax, plm_i, dlm_i)
621 >       call Associated_Legendre(cti, ShapeMap%Shapes(st1)%bigM, &
622 >            ShapeMap%Shapes(st1)%bigL, LMAX, &
623 >            plm_i, dlm_i)
624  
625 <       call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, &
625 >       call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, MMAX, &
626              CHEBYSHEV_TN, tm_i, dtm_i)
627 <       call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, &
627 >       call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, MMAX, &
628              CHEBYSHEV_UN, um_i, dum_i)
629 <      
629 >
630         sigma_i = 0.0d0
631         s_i = 0.0d0
632         eps_i = 0.0d0
# Line 645 | Line 673 | contains  
673               dPhuncdUz = coeff*(spi * dum_i(m-1)*dcpiduz + dspiduz *um_i(m-1))
674            endif
675  
676 <          sigma_i = sigma_i + plm_i(l,m)*Phunc
677 <          
678 <          dsigmaidx = dsigmaidx + plm_i(l,m)*dPhuncdX + &
679 <               Phunc * dlm_i(l,m) * dctidx
680 <          dsigmaidy = dsigmaidy + plm_i(l,m)*dPhuncdY + &
681 <               Phunc * dlm_i(l,m) * dctidy
682 <          dsigmaidz = dsigmaidz + plm_i(l,m)*dPhuncdZ + &
683 <               Phunc * dlm_i(l,m) * dctidz
684 <          
685 <          dsigmaidux = dsigmaidux + plm_i(l,m)* dPhuncdUx + &
686 <               Phunc * dlm_i(l,m) * dctidux
687 <          dsigmaiduy = dsigmaiduy + plm_i(l,m)* dPhuncdUy + &
688 <               Phunc * dlm_i(l,m) * dctiduy
689 <          dsigmaiduz = dsigmaiduz + plm_i(l,m)* dPhuncdUz + &
690 <               Phunc * dlm_i(l,m) * dctiduz
691 <
676 >          sigma_i = sigma_i + plm_i(m,l)*Phunc
677 >          write(*,*) 'dsigmaidux = ', dsigmaidux
678 >          write(*,*) 'Phunc = ', Phunc
679 >          dsigmaidx = dsigmaidx + plm_i(m,l)*dPhuncdX + &
680 >               Phunc * dlm_i(m,l) * dctidx
681 >          dsigmaidy = dsigmaidy + plm_i(m,l)*dPhuncdY + &
682 >               Phunc * dlm_i(m,l) * dctidy
683 >          dsigmaidz = dsigmaidz + plm_i(m,l)*dPhuncdZ + &
684 >               Phunc * dlm_i(m,l) * dctidz
685 >          dsigmaidux = dsigmaidux + plm_i(m,l)* dPhuncdUx + &
686 >               Phunc * dlm_i(m,l) * dctidux
687 >          dsigmaiduy = dsigmaiduy + plm_i(m,l)* dPhuncdUy + &
688 >               Phunc * dlm_i(m,l) * dctiduy
689 >          dsigmaiduz = dsigmaiduz + plm_i(m,l)* dPhuncdUz + &
690 >               Phunc * dlm_i(m,l) * dctiduz
691 >          write(*,*) 'dsigmaidux = ', dsigmaidux, '; dPhuncdUx = ', dPhuncdUx, &
692 >                     '; dctidux = ', dctidux, '; plm_i(m,l) = ', plm_i(m,l), &
693 >                     '; dlm_i(m,l) = ', dlm_i(m,l), '; m = ', m, '; l = ', l
694         end do
695  
696         do lm = 1, ShapeMap%Shapes(st1)%nRangeFuncs
# Line 668 | Line 698 | contains  
698            m = ShapeMap%Shapes(st1)%RangeFuncMValue(lm)
699            coeff = ShapeMap%Shapes(st1)%RangeFuncCoefficient(lm)
700            function_type = ShapeMap%Shapes(st1)%RangeFunctionType(lm)
701 <          
701 >
702            if ((function_type .eq. SH_COS).or.(m.eq.0)) then
703               Phunc = coeff * tm_i(m)
704               dPhuncdX = coeff * dtm_i(m) * dcpidx
# Line 687 | Line 717 | contains  
717               dPhuncdUz = coeff*(spi * dum_i(m-1)*dcpiduz + dspiduz *um_i(m-1))
718            endif
719  
720 <          s_i = s_i + plm_i(l,m)*Phunc
691 <          
692 <          dsidx = dsidx + plm_i(l,m)*dPhuncdX + &
693 <               Phunc * dlm_i(l,m) * dctidx
694 <          dsidy = dsidy + plm_i(l,m)*dPhuncdY + &
695 <               Phunc * dlm_i(l,m) * dctidy
696 <          dsidz = dsidz + plm_i(l,m)*dPhuncdZ + &
697 <               Phunc * dlm_i(l,m) * dctidz
698 <          
699 <          dsidux = dsidux + plm_i(l,m)* dPhuncdUx + &
700 <               Phunc * dlm_i(l,m) * dctidux
701 <          dsiduy = dsiduy + plm_i(l,m)* dPhuncdUy + &
702 <               Phunc * dlm_i(l,m) * dctiduy
703 <          dsiduz = dsiduz + plm_i(l,m)* dPhuncdUz + &
704 <               Phunc * dlm_i(l,m) * dctiduz      
720 >          s_i = s_i + plm_i(m,l)*Phunc
721  
722 +          dsidx = dsidx + plm_i(m,l)*dPhuncdX + &
723 +               Phunc * dlm_i(m,l) * dctidx
724 +          dsidy = dsidy + plm_i(m,l)*dPhuncdY + &
725 +               Phunc * dlm_i(m,l) * dctidy
726 +          dsidz = dsidz + plm_i(m,l)*dPhuncdZ + &
727 +               Phunc * dlm_i(m,l) * dctidz
728 +
729 +          dsidux = dsidux + plm_i(m,l)* dPhuncdUx + &
730 +               Phunc * dlm_i(m,l) * dctidux
731 +          dsiduy = dsiduy + plm_i(m,l)* dPhuncdUy + &
732 +               Phunc * dlm_i(m,l) * dctiduy
733 +          dsiduz = dsiduz + plm_i(m,l)* dPhuncdUz + &
734 +               Phunc * dlm_i(m,l) * dctiduz      
735 +
736         end do
737 <              
737 >
738         do lm = 1, ShapeMap%Shapes(st1)%nStrengthFuncs
739            l = ShapeMap%Shapes(st1)%StrengthFuncLValue(lm)
740            m = ShapeMap%Shapes(st1)%StrengthFuncMValue(lm)
741            coeff = ShapeMap%Shapes(st1)%StrengthFuncCoefficient(lm)
742            function_type = ShapeMap%Shapes(st1)%StrengthFunctionType(lm)
743 <          
743 >
744            if ((function_type .eq. SH_COS).or.(m.eq.0)) then
745               Phunc = coeff * tm_i(m)
746               dPhuncdX = coeff * dtm_i(m) * dcpidx
# Line 729 | Line 759 | contains  
759               dPhuncdUz = coeff*(spi * dum_i(m-1)*dcpiduz + dspiduz *um_i(m-1))
760            endif
761  
762 <          eps_i = eps_i + plm_i(l,m)*Phunc
733 <          
734 <          depsidx = depsidx + plm_i(l,m)*dPhuncdX + &
735 <               Phunc * dlm_i(l,m) * dctidx
736 <          depsidy = depsidy + plm_i(l,m)*dPhuncdY + &
737 <               Phunc * dlm_i(l,m) * dctidy
738 <          depsidz = depsidz + plm_i(l,m)*dPhuncdZ + &
739 <               Phunc * dlm_i(l,m) * dctidz
740 <          
741 <          depsidux = depsidux + plm_i(l,m)* dPhuncdUx + &
742 <               Phunc * dlm_i(l,m) * dctidux
743 <          depsiduy = depsiduy + plm_i(l,m)* dPhuncdUy + &
744 <               Phunc * dlm_i(l,m) * dctiduy
745 <          depsiduz = depsiduz + plm_i(l,m)* dPhuncdUz + &
746 <               Phunc * dlm_i(l,m) * dctiduz      
762 >          eps_i = eps_i + plm_i(m,l)*Phunc
763  
764 +          depsidx = depsidx + plm_i(m,l)*dPhuncdX + &
765 +               Phunc * dlm_i(m,l) * dctidx
766 +          depsidy = depsidy + plm_i(m,l)*dPhuncdY + &
767 +               Phunc * dlm_i(m,l) * dctidy
768 +          depsidz = depsidz + plm_i(m,l)*dPhuncdZ + &
769 +               Phunc * dlm_i(m,l) * dctidz
770 +
771 +          depsidux = depsidux + plm_i(m,l)* dPhuncdUx + &
772 +               Phunc * dlm_i(m,l) * dctidux
773 +          depsiduy = depsiduy + plm_i(m,l)* dPhuncdUy + &
774 +               Phunc * dlm_i(m,l) * dctiduy
775 +          depsiduz = depsiduz + plm_i(m,l)* dPhuncdUz + &
776 +               Phunc * dlm_i(m,l) * dctiduz      
777 +
778         end do
779  
780      endif
751 –      
752 –       ! now do j:
781  
782 +    ! now do j:
783 +
784      if (ShapeMap%Shapes(st2)%isLJ) then
785         sigma_j = ShapeMap%Shapes(st2)%sigma
786         s_j = ShapeMap%Shapes(st2)%sigma
# Line 774 | Line 804 | contains  
804         depsjduy = 0.0d0
805         depsjduz = 0.0d0
806      else
807 <      
807 >
808   #ifdef IS_MPI
809         ! rotate the inter-particle separation into the two different
810         ! body-fixed coordinate systems:
811         ! negative sign because this is the vector from j to i:
812 <      
812 >
813         xj = -(A_Col(1,atom2)*d(1) + A_Col(2,atom2)*d(2) + A_Col(3,atom2)*d(3))
814         yj = -(A_Col(4,atom2)*d(1) + A_Col(5,atom2)*d(2) + A_Col(6,atom2)*d(3))
815         zj = -(A_Col(7,atom2)*d(1) + A_Col(8,atom2)*d(2) + A_Col(9,atom2)*d(3))
# Line 787 | Line 817 | contains  
817         ! rotate the inter-particle separation into the two different
818         ! body-fixed coordinate systems:
819         ! negative sign because this is the vector from j to i:
820 <      
820 >
821         xj = -(a(1,atom2)*d(1) + a(2,atom2)*d(2) + a(3,atom2)*d(3))
822         yj = -(a(4,atom2)*d(1) + a(5,atom2)*d(2) + a(6,atom2)*d(3))
823         zj = -(a(7,atom2)*d(1) + a(8,atom2)*d(2) + a(9,atom2)*d(3))
824   #endif
825 <      
825 >
826         xj2 = xj*xj
827         yj2 = yj*yj
828         zj2 = zj*zj
799 –      
800 –       projj = sqrt(xj2 + yj2)
801 –       projj3 = projj*projj*projj
802 –      
829         ctj = zj / rij
830 +
831 +       if (ctj .gt. 1.0_dp) ctj = 1.0_dp
832 +       if (ctj .lt. -1.0_dp) ctj = -1.0_dp
833 +
834         dctjdx = - zj * xj / r3
835         dctjdy = - zj * yj / r3
836         dctjdz = 1.0d0 / rij - zj2 / r3
837 <       dctjdux =  yj / rij
838 <       dctjduy = -xj / rij
839 <       dctjduz = 0.0d0
840 <      
841 <       cpj = xj / projj
842 <       dcpjdx = 1.0d0 / projj - xj2 / projj3
843 <       dcpjdy = - xj * yj / projj3
844 <       dcpjdz = 0.0d0
845 <       dcpjdux = xj * yj * zj / projj3
846 <       dcpjduy = -zj * (1.0d0 / projj - xj2 / projj3)
847 <       dcpjduz = -yj * (1.0d0 / projj - xj2 / projj3)  - (xj2 * yj / projj3)
848 <      
849 <       spj = yj / projj
850 <       dspjdx = - xj * yj / projj3
851 <       dspjdy = 1.0d0 / projj - yj2 / projj3
837 >       dctjdux = - (zi * xj2) / r3
838 >       dctjduy = - (zj * yj2) / r3
839 >       dctjduz = zj / rij - (zj2 * zj) / r3
840 >
841 >       ! this is an attempt to try to truncate the singularity when
842 >       ! sin(theta) is near 0.0:
843 >
844 >       stj2 = 1.0_dp - ctj*ctj
845 >       if (dabs(stj2) .lt. 1.0d-12) then
846 >          projj = sqrt(rij * 1.0d-12)
847 >          dcpjdx = 1.0d0 / projj
848 >          dcpjdy = 0.0d0
849 >          dcpjdux = xj / projj
850 >          dcpjduy = 0.0d0
851 >          dspjdx = 0.0d0
852 >          dspjdy = 1.0d0 / projj
853 >          dspjdux = 0.0d0
854 >          dspjduy = yj / projj
855 >       else
856 >          projj = sqrt(xj2 + yj2)
857 >          projj3 = projj*projj*projj
858 >          dcpjdx = 1.0d0 / projj - xj2 / projj3
859 >          dcpjdy = - xj * yj / projj3
860 >          dcpjdux = xj / projj - (xj2 * xj) / projj3
861 >          dcpjduy = - (xj * yj2) / projj3
862 >          dspjdx = - xj * yj / projj3
863 >          dspjdy = 1.0d0 / projj - yj2 / projj3
864 >          dspjdux = - (yj * xj2) / projj3
865 >          dspjduy = yj / projj - (yj2 * yj) / projj3
866 >       endif
867 >
868 >       cpj = xj / projj
869 >       dcpjdz = 0.0d0
870 >       dcpjduz = 0.0d0
871 >
872 >       spj = yj / projj
873         dspjdz = 0.0d0
874 <       dspjdux = -zj * (1.0d0 / projj - yj2 / projj3)
875 <       dspjduy = xj * yj * zj / projj3
876 <       dspjduz = xj * (1.0d0 / projj - yi2 / projj3) + (xj * yj2 / projj3)
877 <      
878 <       call Associated_Legendre(ctj, ShapeMap%Shapes(st2)%bigL, &
879 <            ShapeMap%Shapes(st2)%bigM, lmax, plm_j, dlm_j)
880 <      
881 <       call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, &
874 >       dspjduz = 0.0d0
875 >
876 >
877 > !       write(*,*) 'dcpdu = ' ,dcpidux, dcpiduy, dcpiduz
878 > !       write(*,*) 'dcpdu = ' ,dcpjdux, dcpjduy, dcpjduz
879 >       call Associated_Legendre(ctj, ShapeMap%Shapes(st2)%bigM, &
880 >            ShapeMap%Shapes(st2)%bigL, LMAX, &
881 >            plm_j, dlm_j)
882 >
883 >       call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, MMAX, &
884              CHEBYSHEV_TN, tm_j, dtm_j)
885 <       call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, &
885 >       call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, MMAX, &
886              CHEBYSHEV_UN, um_j, dum_j)
887 <      
887 >
888         sigma_j = 0.0d0
889         s_j = 0.0d0
890         eps_j = 0.0d0
# Line 877 | Line 930 | contains  
930               dPhuncdUy = coeff*(spj * dum_j(m-1)*dcpjduy + dspjduy *um_j(m-1))
931               dPhuncdUz = coeff*(spj * dum_j(m-1)*dcpjduz + dspjduz *um_j(m-1))
932            endif
880 –
881 –          sigma_j = sigma_j + plm_j(l,m)*Phunc
882 –          
883 –          dsigmajdx = dsigmajdx + plm_j(l,m)*dPhuncdX + &
884 –               Phunc * dlm_j(l,m) * dctjdx
885 –          dsigmajdy = dsigmajdy + plm_j(l,m)*dPhuncdY + &
886 –               Phunc * dlm_j(l,m) * dctjdy
887 –          dsigmajdz = dsigmajdz + plm_j(l,m)*dPhuncdZ + &
888 –               Phunc * dlm_j(l,m) * dctjdz
889 –          
890 –          dsigmajdux = dsigmajdux + plm_j(l,m)* dPhuncdUx + &
891 –               Phunc * dlm_j(l,m) * dctjdux
892 –          dsigmajduy = dsigmajduy + plm_j(l,m)* dPhuncdUy + &
893 –               Phunc * dlm_j(l,m) * dctjduy
894 –          dsigmajduz = dsigmajduz + plm_j(l,m)* dPhuncdUz + &
895 –               Phunc * dlm_j(l,m) * dctjduz
933  
934 +          sigma_j = sigma_j + plm_j(m,l)*Phunc
935 +
936 +          dsigmajdx = dsigmajdx + plm_j(m,l)*dPhuncdX + &
937 +               Phunc * dlm_j(m,l) * dctjdx
938 +          dsigmajdy = dsigmajdy + plm_j(m,l)*dPhuncdY + &
939 +               Phunc * dlm_j(m,l) * dctjdy
940 +          dsigmajdz = dsigmajdz + plm_j(m,l)*dPhuncdZ + &
941 +               Phunc * dlm_j(m,l) * dctjdz
942 +
943 +          dsigmajdux = dsigmajdux + plm_j(m,l)* dPhuncdUx + &
944 +               Phunc * dlm_j(m,l) * dctjdux
945 +          dsigmajduy = dsigmajduy + plm_j(m,l)* dPhuncdUy + &
946 +               Phunc * dlm_j(m,l) * dctjduy
947 +          dsigmajduz = dsigmajduz + plm_j(m,l)* dPhuncdUz + &
948 +               Phunc * dlm_j(m,l) * dctjduz
949 +
950         end do
951  
952         do lm = 1, ShapeMap%Shapes(st2)%nRangeFuncs
# Line 920 | Line 973 | contains  
973               dPhuncdUz = coeff*(spj * dum_j(m-1)*dcpjduz + dspjduz *um_j(m-1))
974            endif
975  
976 <          s_j = s_j + plm_j(l,m)*Phunc
977 <          
978 <          dsjdx = dsjdx + plm_j(l,m)*dPhuncdX + &
979 <               Phunc * dlm_j(l,m) * dctjdx
980 <          dsjdy = dsjdy + plm_j(l,m)*dPhuncdY + &
981 <               Phunc * dlm_j(l,m) * dctjdy
982 <          dsjdz = dsjdz + plm_j(l,m)*dPhuncdZ + &
983 <               Phunc * dlm_j(l,m) * dctjdz
984 <          
985 <          dsjdux = dsjdux + plm_j(l,m)* dPhuncdUx + &
986 <               Phunc * dlm_j(l,m) * dctjdux
987 <          dsjduy = dsjduy + plm_j(l,m)* dPhuncdUy + &
988 <               Phunc * dlm_j(l,m) * dctjduy
989 <          dsjduz = dsjduz + plm_j(l,m)* dPhuncdUz + &
990 <               Phunc * dlm_j(l,m) * dctjduz
976 >          s_j = s_j + plm_j(m,l)*Phunc
977 >
978 >          dsjdx = dsjdx + plm_j(m,l)*dPhuncdX + &
979 >               Phunc * dlm_j(m,l) * dctjdx
980 >          dsjdy = dsjdy + plm_j(m,l)*dPhuncdY + &
981 >               Phunc * dlm_j(m,l) * dctjdy
982 >          dsjdz = dsjdz + plm_j(m,l)*dPhuncdZ + &
983 >               Phunc * dlm_j(m,l) * dctjdz
984 >
985 >          dsjdux = dsjdux + plm_j(m,l)* dPhuncdUx + &
986 >               Phunc * dlm_j(m,l) * dctjdux
987 >          dsjduy = dsjduy + plm_j(m,l)* dPhuncdUy + &
988 >               Phunc * dlm_j(m,l) * dctjduy
989 >          dsjduz = dsjduz + plm_j(m,l)* dPhuncdUz + &
990 >               Phunc * dlm_j(m,l) * dctjduz
991  
992         end do
993  
# Line 962 | Line 1015 | contains  
1015               dPhuncdUz = coeff*(spj * dum_j(m-1)*dcpjduz + dspjduz *um_j(m-1))
1016            endif
1017  
1018 <          eps_j = eps_j + plm_j(l,m)*Phunc
966 <          
967 <          depsjdx = depsjdx + plm_j(l,m)*dPhuncdX + &
968 <               Phunc * dlm_j(l,m) * dctjdx
969 <          depsjdy = depsjdy + plm_j(l,m)*dPhuncdY + &
970 <               Phunc * dlm_j(l,m) * dctjdy
971 <          depsjdz = depsjdz + plm_j(l,m)*dPhuncdZ + &
972 <               Phunc * dlm_j(l,m) * dctjdz
973 <          
974 <          depsjdux = depsjdux + plm_j(l,m)* dPhuncdUx + &
975 <               Phunc * dlm_j(l,m) * dctjdux
976 <          depsjduy = depsjduy + plm_j(l,m)* dPhuncdUy + &
977 <               Phunc * dlm_j(l,m) * dctjduy
978 <          depsjduz = depsjduz + plm_j(l,m)* dPhuncdUz + &
979 <               Phunc * dlm_j(l,m) * dctjduz
1018 > !          write(*,*) 'l,m = ', l, m, coeff, dPhuncdUx, dPhuncdUy, dPhuncdUz
1019  
1020 +          eps_j = eps_j + plm_j(m,l)*Phunc
1021 +
1022 +          depsjdx = depsjdx + plm_j(m,l)*dPhuncdX + &
1023 +               Phunc * dlm_j(m,l) * dctjdx
1024 +          depsjdy = depsjdy + plm_j(m,l)*dPhuncdY + &
1025 +               Phunc * dlm_j(m,l) * dctjdy
1026 +          depsjdz = depsjdz + plm_j(m,l)*dPhuncdZ + &
1027 +               Phunc * dlm_j(m,l) * dctjdz
1028 +
1029 +          depsjdux = depsjdux + plm_j(m,l)* dPhuncdUx + &
1030 +               Phunc * dlm_j(m,l) * dctjdux
1031 +          depsjduy = depsjduy + plm_j(m,l)* dPhuncdUy + &
1032 +               Phunc * dlm_j(m,l) * dctjduy
1033 +          depsjduz = depsjduz + plm_j(m,l)* dPhuncdUz + &
1034 +               Phunc * dlm_j(m,l) * dctjduz
1035 +
1036         end do
1037  
1038      endif
# Line 985 | Line 1040 | contains  
1040      ! phew, now let's assemble the potential energy:
1041  
1042      sigma = 0.5*(sigma_i + sigma_j)
1043 <
1043 > !    write(*,*) sigma_i, ' = sigma_i; ', sigma_j, ' = sigma_j'
1044      dsigmadxi = 0.5*dsigmaidx
1045      dsigmadyi = 0.5*dsigmaidy
1046      dsigmadzi = 0.5*dsigmaidz
# Line 1017 | Line 1072 | contains  
1072      dsduzj = 0.5*dsjduz
1073  
1074      eps = sqrt(eps_i * eps_j)
1075 <
1075 >    write(*,*) 'dsidu = ', dsidux, dsiduy, dsiduz
1076 >    write(*,*) 'dsigidu = ', dsigmaidux, dsigmaiduy, dsigmaiduz
1077 > !    write(*,*) sigma_i, ' is sigma i; ', s_i, ' is s i; ', eps_i, ' is eps i'
1078      depsdxi = eps_j * depsidx / (2.0d0 * eps)
1079      depsdyi = eps_j * depsidy / (2.0d0 * eps)
1080      depsdzi = eps_j * depsidz / (2.0d0 * eps)
# Line 1031 | Line 1088 | contains  
1088      depsduxj = eps_i * depsjdux / (2.0d0 * eps)
1089      depsduyj = eps_i * depsjduy / (2.0d0 * eps)
1090      depsduzj = eps_i * depsjduz / (2.0d0 * eps)
1091 <    
1091 >
1092 > !    write(*,*) 'depsidu = ', depsidux, depsiduy, depsiduz
1093 > !    write(*,*) 'depsjdu = ', depsjdux, depsjduy, depsjduz
1094 >
1095 > !    write(*,*) 'depsdui = ', depsduxi, depsduyi, depsduzi
1096 > !    write(*,*) 'depsduj = ', depsduxj, depsduyj, depsduzj
1097 > !!$
1098 > !!$    write(*,*) 's, sig, eps = ', s, sigma, eps
1099 >
1100      rtdenom = rij-sigma+s
1101      rt = s / rtdenom
1102  
# Line 1039 | Line 1104 | contains  
1104      drtdyi = (dsdyi + rt * (drdyi - dsigmadyi + dsdyi)) / rtdenom
1105      drtdzi = (dsdzi + rt * (drdzi - dsigmadzi + dsdzi)) / rtdenom
1106      drtduxi = (dsduxi + rt * (drduxi - dsigmaduxi + dsduxi)) / rtdenom
1107 +    write(*,*) dsduxi, ' is dsduxi; ', drduxi, ' is drduxi; ', dsigmaduxi, &
1108 +               ' is dsigmaduxi; ', dsduxi, ' is dsduxi'
1109      drtduyi = (dsduyi + rt * (drduyi - dsigmaduyi + dsduyi)) / rtdenom
1110      drtduzi = (dsduzi + rt * (drduzi - dsigmaduzi + dsduzi)) / rtdenom
1111      drtdxj = (dsdxj + rt * (drdxj - dsigmadxj + dsdxj)) / rtdenom
# Line 1047 | Line 1114 | contains  
1114      drtduxj = (dsduxj + rt * (drduxj - dsigmaduxj + dsduxj)) / rtdenom
1115      drtduyj = (dsduyj + rt * (drduyj - dsigmaduyj + dsduyj)) / rtdenom
1116      drtduzj = (dsduzj + rt * (drduzj - dsigmaduzj + dsduzj)) / rtdenom
1117 <    
1117 >
1118 > !    write(*,*) 'drtd_i = ', drtdxi, drtdyi, drtdzi
1119 > !    write(*,*) 'drtdu_i = ', drtduxi, drtduyi, drtduzi
1120 >
1121      rt2 = rt*rt
1122      rt3 = rt2*rt
1123      rt5 = rt2*rt3
# Line 1056 | Line 1126 | contains  
1126      rt12 = rt6*rt6
1127      rt126 = rt12 - rt6
1128  
1129 +    pot_temp = 4.0d0 * eps * rt126
1130 +
1131 +    vpair = vpair + pot_temp
1132      if (do_pot) then
1133   #ifdef IS_MPI
1134 <       pot_row(atom1) = pot_row(atom1) + 2.0d0*eps*rt126*sw
1135 <       pot_col(atom2) = pot_col(atom2) + 2.0d0*eps*rt126*sw
1134 >       pot_row(atom1) = pot_row(atom1) + 0.5d0*pot_temp*sw
1135 >       pot_col(atom2) = pot_col(atom2) + 0.5d0*pot_temp*sw
1136   #else
1137 <       pot = pot + 4.0d0*eps*rt126*sw
1137 >       pot = pot + pot_temp*sw
1138   #endif
1139      endif
1140 <    
1140 >
1141 > !!$    write(*,*) 'drtdu, depsdu = ', drtduxi, depsduxi
1142 >
1143      dvdxi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdxi + 4.0d0*depsdxi*rt126
1144      dvdyi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdyi + 4.0d0*depsdyi*rt126
1145      dvdzi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdzi + 4.0d0*depsdzi*rt126
# Line 1078 | Line 1153 | contains  
1153      dvduxj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduxj + 4.0d0*depsduxj*rt126
1154      dvduyj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduyj + 4.0d0*depsduyj*rt126
1155      dvduzj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduzj + 4.0d0*depsduzj*rt126
1156 <
1156 >    write(*,*) 'drtduxi = ', drtduxi, ' depsduxi = ', depsduxi
1157      ! do the torques first since they are easy:
1158      ! remember that these are still in the body fixed axes
1159  
1160 +
1161 + !!$    write(*,*) 'sw = ', sw
1162 + !!$    write(*,*) 'dvdu1 = ', dvduxi, dvduyi, dvduzi
1163 + !!$    write(*,*) 'dvdu2 = ', dvduxj, dvduyj, dvduzj
1164 + !!$
1165 + !    txi =  (dvduzi - dvduyi) * sw
1166 + !    tyi =  (dvduxi - dvduzi) * sw
1167 + !    tzi =  (dvduyi - dvduxi) * sw
1168 +
1169 + !    txj = (dvduzj - dvduyj) * sw
1170 + !    tyj = (dvduxj - dvduzj) * sw
1171 + !    tzj = (dvduyj - dvduxj) * sw
1172 +
1173      txi = dvduxi * sw
1174      tyi = dvduyi * sw
1175      tzi = dvduzi * sw
# Line 1090 | Line 1178 | contains  
1178      tyj = dvduyj * sw
1179      tzj = dvduzj * sw
1180  
1181 +    write(*,*) 't1 = ', txi, tyi, tzi
1182 +    write(*,*) 't2 = ', txj, tyj, tzj
1183 +
1184      ! go back to lab frame using transpose of rotation matrix:
1185 <    
1185 >
1186   #ifdef IS_MPI
1187      t_Row(1,atom1) = t_Row(1,atom1) + a_Row(1,atom1)*txi + &
1188           a_Row(4,atom1)*tyi + a_Row(7,atom1)*tzi
# Line 1099 | Line 1190 | contains  
1190           a_Row(5,atom1)*tyi + a_Row(8,atom1)*tzi
1191      t_Row(3,atom1) = t_Row(3,atom1) + a_Row(3,atom1)*txi + &
1192           a_Row(6,atom1)*tyi + a_Row(9,atom1)*tzi
1193 <    
1193 >
1194      t_Col(1,atom2) = t_Col(1,atom2) + a_Col(1,atom2)*txj + &
1195           a_Col(4,atom2)*tyj + a_Col(7,atom2)*tzj
1196      t_Col(2,atom2) = t_Col(2,atom2) + a_Col(2,atom2)*txj + &
1197 <            a_Col(5,atom2)*tyj + a_Col(8,atom2)*tzj
1197 >         a_Col(5,atom2)*tyj + a_Col(8,atom2)*tzj
1198      t_Col(3,atom2) = t_Col(3,atom2) + a_Col(3,atom2)*txj + &
1199           a_Col(6,atom2)*tyj + a_Col(9,atom2)*tzj
1200   #else
1201      t(1,atom1) = t(1,atom1) + a(1,atom1)*txi + a(4,atom1)*tyi + a(7,atom1)*tzi
1202      t(2,atom1) = t(2,atom1) + a(2,atom1)*txi + a(5,atom1)*tyi + a(8,atom1)*tzi
1203      t(3,atom1) = t(3,atom1) + a(3,atom1)*txi + a(6,atom1)*tyi + a(9,atom1)*tzi
1204 <    
1204 >
1205      t(1,atom2) = t(1,atom2) + a(1,atom2)*txj + a(4,atom2)*tyj + a(7,atom2)*tzj
1206      t(2,atom2) = t(2,atom2) + a(2,atom2)*txj + a(5,atom2)*tyj + a(8,atom2)*tzj
1207      t(3,atom2) = t(3,atom2) + a(3,atom2)*txj + a(6,atom2)*tyj + a(9,atom2)*tzj
1208   #endif
1209      ! Now, on to the forces:
1210 <    
1210 >
1211      ! first rotate the i terms back into the lab frame:
1212 <    
1212 >
1213      fxi = dvdxi * sw
1214      fyi = dvdyi * sw
1215      fzi = dvdzi * sw
# Line 1127 | Line 1218 | contains  
1218      fyj = dvdyj * sw
1219      fzj = dvdzj * sw
1220  
1221 +
1222   #ifdef IS_MPI
1223      fxii = a_Row(1,atom1)*fxi + a_Row(4,atom1)*fyi + a_Row(7,atom1)*fzi
1224      fyii = a_Row(2,atom1)*fxi + a_Row(5,atom1)*fyi + a_Row(8,atom1)*fzi
# Line 1139 | Line 1231 | contains  
1231      fxii = a(1,atom1)*fxi + a(4,atom1)*fyi + a(7,atom1)*fzi
1232      fyii = a(2,atom1)*fxi + a(5,atom1)*fyi + a(8,atom1)*fzi
1233      fzii = a(3,atom1)*fxi + a(6,atom1)*fyi + a(9,atom1)*fzi
1234 <    
1234 >
1235      fxjj = a(1,atom2)*fxj + a(4,atom2)*fyj + a(7,atom2)*fzj
1236      fyjj = a(2,atom2)*fxj + a(5,atom2)*fyj + a(8,atom2)*fzj
1237      fzjj = a(3,atom2)*fxj + a(6,atom2)*fyj + a(9,atom2)*fzj
# Line 1148 | Line 1240 | contains  
1240      fxij = -fxii
1241      fyij = -fyii
1242      fzij = -fzii
1243 <    
1243 >
1244      fxji = -fxjj
1245      fyji = -fyjj
1246      fzji = -fzjj
1247  
1248 <    fxradial = fxii + fxji
1249 <    fyradial = fyii + fyji
1250 <    fzradial = fzii + fzji
1251 <
1248 >    fxradial = 0.5_dp * (fxii + fxji)
1249 >    fyradial = 0.5_dp * (fyii + fyji)
1250 >    fzradial = 0.5_dp * (fzii + fzji)
1251 >    write(*,*) fxradial, ' is fxrad; ', fyradial, ' is fyrad; ', fzradial, 'is fzrad'
1252   #ifdef IS_MPI
1253      f_Row(1,atom1) = f_Row(1,atom1) + fxradial
1254      f_Row(2,atom1) = f_Row(2,atom1) + fyradial
1255      f_Row(3,atom1) = f_Row(3,atom1) + fzradial
1256 <    
1256 >
1257      f_Col(1,atom2) = f_Col(1,atom2) - fxradial
1258      f_Col(2,atom2) = f_Col(2,atom2) - fyradial
1259      f_Col(3,atom2) = f_Col(3,atom2) - fzradial
# Line 1169 | Line 1261 | contains  
1261      f(1,atom1) = f(1,atom1) + fxradial
1262      f(2,atom1) = f(2,atom1) + fyradial
1263      f(3,atom1) = f(3,atom1) + fzradial
1264 <    
1264 >
1265      f(1,atom2) = f(1,atom2) - fxradial
1266      f(2,atom2) = f(2,atom2) - fyradial
1267      f(3,atom2) = f(3,atom2) - fzradial
# Line 1182 | Line 1274 | contains  
1274      id1 = atom1
1275      id2 = atom2
1276   #endif
1277 <    
1277 >
1278      if (molMembershipList(id1) .ne. molMembershipList(id2)) then
1279 <      
1279 >
1280         fpair(1) = fpair(1) + fxradial
1281         fpair(2) = fpair(2) + fyradial
1282         fpair(3) = fpair(3) + fzradial
1283 <      
1283 >
1284      endif
1285 <    
1285 >
1286    end subroutine do_shape_pair
1287 <    
1288 <  SUBROUTINE Associated_Legendre(x, l, m, lmax, plm, dlm)
1289 <    
1287 >
1288 >  SUBROUTINE Associated_Legendre(x, l, m, lmax, plm, dlm)        
1289 >
1290      ! Purpose: Compute the associated Legendre functions
1291      !          Plm(x) and their derivatives Plm'(x)
1292      ! Input :  x  --- Argument of Plm(x)
# Line 1210 | Line 1302 | contains  
1302      !
1303      ! The original Fortran77 codes can be found here:
1304      ! http://iris-lee3.ece.uiuc.edu/~jjin/routines/routines.html
1305 <    
1306 <    real (kind=8), intent(in) :: x
1305 >
1306 >    real (kind=dp), intent(in) :: x
1307      integer, intent(in) :: l, m, lmax
1308 <    real (kind=8), dimension(0:lmax,0:m), intent(out) :: PLM, DLM
1308 >    real (kind=dp), dimension(0:lmax,0:m), intent(out) :: PLM, DLM
1309      integer :: i, j, ls
1310 <    real (kind=8) :: xq, xs
1310 >    real (kind=dp) :: xq, xs
1311  
1312      ! zero out both arrays:
1313      DO I = 0, m
1314         DO J = 0, l
1315 <          PLM(J,I) = 0.0D0
1316 <          DLM(J,I) = 0.0D0
1315 >          PLM(J,I) = 0.0_dp
1316 >          DLM(J,I) = 0.0_dp
1317         end DO
1318      end DO
1319  
1320      ! start with 0,0:
1321      PLM(0,0) = 1.0D0
1322 <  
1322 >
1323      ! x = +/- 1 functions are easy:
1324      IF (abs(X).EQ.1.0D0) THEN
1325         DO I = 1, m
# Line 1254 | Line 1346 | contains  
1346      DO I = 1, l
1347         PLM(I, I) = -LS*(2.0D0*I-1.0D0)*XQ*PLM(I-1, I-1)
1348      enddo
1349 <    
1349 >
1350      DO I = 0, l
1351         PLM(I, I+1)=(2.0D0*I+1.0D0)*X*PLM(I, I)
1352      enddo
1353 <    
1353 >
1354      DO I = 0, l
1355         DO J = I+2, m
1356            PLM(I, J)=((2.0D0*J-1.0D0)*X*PLM(I,J-1) - &
1357                 (I+J-1.0D0)*PLM(I,J-2))/(J-I)
1358         end DO
1359      end DO
1360 <  
1360 >
1361      DLM(0, 0)=0.0D0
1270 –    
1362      DO J = 1, m
1363         DLM(0, J)=LS*J*(PLM(0,J-1)-X*PLM(0,J))/XS
1364      end DO
1365 <    
1365 >
1366      DO I = 1, l
1367         DO J = I, m
1368            DLM(I,J) = LS*I*X*PLM(I, J)/XS + (J+I)*(J-I+1.0D0)/XQ*PLM(I-1, J)
1369         end DO
1370      end DO
1371 <    
1371 >
1372      RETURN
1373    END SUBROUTINE Associated_Legendre
1374  
1375  
1376 <  subroutine Orthogonal_Polynomial(x, m, function_type, pl, dpl)
1377 <  
1376 >  subroutine Orthogonal_Polynomial(x, m, mmax, function_type, pl, dpl)
1377 >
1378      ! Purpose: Compute orthogonal polynomials: Tn(x) or Un(x),
1379      !          or Ln(x) or Hn(x), and their derivatives
1380      ! Input :  function_type --- Function code
# Line 1302 | Line 1393 | contains  
1393      !
1394      ! The original Fortran77 codes can be found here:
1395      ! http://iris-lee3.ece.uiuc.edu/~jjin/routines/routines.html
1396 <  
1396 >
1397      real(kind=8), intent(in) :: x
1398 <    integer, intent(in):: m
1398 >    integer, intent(in):: m, mmax
1399      integer, intent(in):: function_type
1400 <    real(kind=8), dimension(0:m), intent(inout) :: pl, dpl
1401 <  
1400 >    real(kind=8), dimension(0:mmax), intent(inout) :: pl, dpl
1401 >
1402      real(kind=8) :: a, b, c, y0, y1, dy0, dy1, yn, dyn
1403      integer :: k
1404  
# Line 1350 | Line 1441 | contains  
1441         DY0 = DY1
1442         DY1 = DYN
1443      end DO
1444 +
1445 +
1446      RETURN
1447 <    
1447 >
1448    end subroutine Orthogonal_Polynomial
1356 –  
1357 – end module shapes
1449  
1450 < subroutine makeShape(nContactFuncs, ContactFuncLValue, &
1451 <     ContactFuncMValue, ContactFunctionType, ContactFuncCoefficient, &
1361 <     nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
1362 <     RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
1363 <     StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
1364 <     myAtid, status)
1450 >  subroutine deallocateShapes(this)
1451 >    type(Shape), pointer :: this
1452  
1453 <  use definitions
1454 <  use shapes, only: newShapeType
1455 <  
1456 <  integer :: nContactFuncs
1370 <  integer :: nRangeFuncs
1371 <  integer :: nStrengthFuncs
1372 <  integer :: status
1373 <  integer :: myAtid
1374 <  
1375 <  integer, dimension(nContactFuncs) :: ContactFuncLValue          
1376 <  integer, dimension(nContactFuncs) :: ContactFuncMValue          
1377 <  integer, dimension(nContactFuncs) :: ContactFunctionType        
1378 <  real(kind=dp), dimension(nContactFuncs) :: ContactFuncCoefficient
1379 <  integer, dimension(nRangeFuncs) :: RangeFuncLValue            
1380 <  integer, dimension(nRangeFuncs) :: RangeFuncMValue            
1381 <  integer, dimension(nRangeFuncs) :: RangeFunctionType          
1382 <  real(kind=dp), dimension(nRangeFuncs) :: RangeFuncCoefficient  
1383 <  integer, dimension(nStrengthFuncs) :: StrengthFuncLValue          
1384 <  integer, dimension(nStrengthFuncs) :: StrengthFuncMValue          
1385 <  integer, dimension(nStrengthFuncs) :: StrengthFunctionType        
1386 <  real(kind=dp), dimension(nStrengthFuncs) :: StrengthFuncCoefficient
1387 <  
1388 <  call newShapeType(nContactFuncs, ContactFuncLValue, &
1389 <       ContactFuncMValue, ContactFunctionType, ContactFuncCoefficient, &
1390 <       nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
1391 <       RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
1392 <       StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
1393 <       myAtid, status)
1453 >    if (associated( this%ContactFuncLValue)) then
1454 >       deallocate(this%ContactFuncLValue)
1455 >       this%ContactFuncLValue => null()
1456 >    end if
1457  
1458 <  return
1459 < end subroutine makeShape
1458 >    if (associated( this%ContactFuncMValue)) then
1459 >       deallocate( this%ContactFuncMValue)
1460 >       this%ContactFuncMValue => null()
1461 >    end if
1462 >    if (associated( this%ContactFunctionType)) then
1463 >       deallocate(this%ContactFunctionType)
1464 >       this%ContactFunctionType => null()
1465 >    end if
1466 >
1467 >    if (associated( this%ContactFuncCoefficient)) then
1468 >       deallocate(this%ContactFuncCoefficient)
1469 >       this%ContactFuncCoefficient => null()
1470 >    end if
1471 >
1472 >    if (associated( this%RangeFuncLValue)) then
1473 >       deallocate(this%RangeFuncLValue)
1474 >       this%RangeFuncLValue => null()
1475 >    end if
1476 >    if (associated( this%RangeFuncMValue)) then
1477 >       deallocate( this%RangeFuncMValue)
1478 >       this%RangeFuncMValue => null()
1479 >    end if
1480 >
1481 >    if (associated( this%RangeFunctionType)) then
1482 >       deallocate( this%RangeFunctionType)
1483 >       this%RangeFunctionType => null()
1484 >    end if
1485 >    if (associated( this%RangeFuncCoefficient)) then
1486 >       deallocate(this%RangeFuncCoefficient)
1487 >       this%RangeFuncCoefficient => null()
1488 >    end if
1489 >
1490 >    if (associated( this%StrengthFuncLValue)) then
1491 >       deallocate(this%StrengthFuncLValue)
1492 >       this%StrengthFuncLValue => null()
1493 >    end if
1494 >
1495 >    if (associated( this%StrengthFuncMValue )) then
1496 >       deallocate(this%StrengthFuncMValue)
1497 >       this%StrengthFuncMValue => null()
1498 >    end if
1499 >
1500 >    if(associated( this%StrengthFunctionType)) then
1501 >       deallocate(this%StrengthFunctionType)
1502 >       this%StrengthFunctionType => null()
1503 >    end if
1504 >    if (associated( this%StrengthFuncCoefficient )) then
1505 >       deallocate(this%StrengthFuncCoefficient)
1506 >       this%StrengthFuncCoefficient => null()
1507 >    end if
1508 >  end subroutine deallocateShapes
1509 >
1510 >  subroutine destroyShapeTypes
1511 >    integer :: i
1512 >    type(Shape), pointer :: thisShape
1513 >
1514 >    ! First walk through and kill the shape
1515 >    do i = 1,ShapeMap%n_shapes
1516 >       thisShape => ShapeMap%Shapes(i)
1517 >       call deallocateShapes(thisShape)
1518 >    end do
1519 >
1520 >    ! set shape map to starting values
1521 >    ShapeMap%n_shapes = 0
1522 >    ShapeMap%currentShape = 0
1523 >
1524 >    if (associated(ShapeMap%Shapes)) then
1525 >       deallocate(ShapeMap%Shapes)
1526 >       ShapeMap%Shapes => null()
1527 >    end if
1528 >
1529 >    if (associated(ShapeMap%atidToShape)) then
1530 >       deallocate(ShapeMap%atidToShape)
1531 >       ShapeMap%atidToShape => null()
1532 >    end if
1533 >
1534 >
1535 >  end subroutine destroyShapeTypes
1536 >
1537 >
1538 > end module shapes

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines