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 2204 by gezelter, Fri Apr 15 22:04:00 2005 UTC vs.
Revision 2756 by gezelter, Wed May 17 15:37:15 2006 UTC

# Line 55 | Line 55 | module shapes
55    implicit none
56  
57    PRIVATE
58 <
58 > #define __FORTRAN90
59 > #include "UseTheForce/DarkSide/fInteractionMap.h"
60    INTEGER, PARAMETER:: CHEBYSHEV_TN = 1
61    INTEGER, PARAMETER:: CHEBYSHEV_UN = 2
62    INTEGER, PARAMETER:: LAGUERRE     = 3
# Line 69 | Line 70 | module shapes
70    public :: newShapeType
71    public :: complete_Shape_FF
72    public :: destroyShapeTypes
73 +  public :: getShapeCut
74  
75    type, private :: Shape
76       integer :: atid
# Line 97 | Line 99 | module shapes
99    type, private :: ShapeList
100       integer :: n_shapes = 0
101       integer :: currentShape = 0
102 <     type (Shape), pointer :: Shapes(:)      => null()
103 <     integer, pointer      :: atidToShape(:) => null()
102 >     type(Shape), pointer :: Shapes(:)      => null()
103 >     integer, pointer     :: atidToShape(:) => null()
104    end type ShapeList
105 <
105 >  
106    type(ShapeList), save :: ShapeMap
107 <
107 >  
108    integer :: lmax
109 <
109 >  
110   contains  
111 <
111 >  
112    subroutine newShapeType(nContactFuncs, ContactFuncLValue, &
113         ContactFuncMValue, ContactFunctionType, ContactFuncCoefficient, &
114         nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
115         RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
116         StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
117 <       myATID, status)
118 <
117 >       c_ident, status)
118 >    
119      integer :: nContactFuncs
120      integer :: nRangeFuncs
121      integer :: nStrengthFuncs
122      integer :: shape_ident
123      integer :: status
124 +    integer :: c_ident
125      integer :: myATID
126      integer :: bigL
127      integer :: bigM
# Line 153 | Line 156 | contains  
156         allocate(ShapeMap%Shapes(nShapeTypes + nLJTypes))
157  
158         ntypes = getSize(atypes)
159 <
160 <       allocate(ShapeMap%atidToShape(0:ntypes))
159 >
160 >       allocate(ShapeMap%atidToShape(ntypes))
161      end if
162  
163      ShapeMap%currentShape = ShapeMap%currentShape + 1
# Line 167 | Line 170 | contains  
170         return
171      endif
172  
173 <    call getElementProperty(atypes, myATID, 'c_ident', me)
173 >    myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
174  
175 <    ShapeMap%atidToShape(me)                         = current
176 <    ShapeMap%Shapes(current)%atid                    = me
175 >    ShapeMap%atidToShape(myATID)                     = current
176 >    ShapeMap%Shapes(current)%atid                    = myATID
177      ShapeMap%Shapes(current)%nContactFuncs           = nContactFuncs
178      ShapeMap%Shapes(current)%nRangeFuncs             = nRangeFuncs
179      ShapeMap%Shapes(current)%nStrengthFuncs          = nStrengthFuncs
# Line 336 | Line 339 | contains  
339      integer :: status
340      integer :: i, j, l, m, lm, function_type
341      real(kind=dp) :: thisDP, sigma
342 <    integer :: alloc_stat, iTheta, iPhi, nSteps, nAtypes, thisIP, current
342 >    integer :: alloc_stat, iTheta, iPhi, nSteps, nAtypes, myATID, current
343      logical :: thisProperty
344  
345      status = 0
# Line 351 | Line 354 | contains  
354      if (nAtypes == 0) then
355         status = -1
356         return
357 <    end if
357 >    end if      
358  
359      ! atypes comes from c side
360 <    do i = 0, nAtypes
361 <
362 <       call getElementProperty(atypes, i, "is_LennardJones", thisProperty)
363 <
360 >    do i = 1, nAtypes
361 >    
362 >       myATID = getFirstMatchingElement(atypes, 'c_ident', i)
363 >       call getElementProperty(atypes, myATID, "is_LennardJones", thisProperty)
364 >        
365         if (thisProperty) then
362 –
366            ShapeMap%currentShape = ShapeMap%currentShape + 1
367            current = ShapeMap%currentShape
368  
369 <          call getElementProperty(atypes, i, "c_ident",  thisIP)
370 <          ShapeMap%atidToShape(thisIP) = current
368 <          ShapeMap%Shapes(current)%atid = thisIP
369 >          ShapeMap%atidToShape(myATID) = current
370 >          ShapeMap%Shapes(current)%atid = myATID
371  
372            ShapeMap%Shapes(current)%isLJ = .true.
373  
374 <          ShapeMap%Shapes(current)%epsilon = getEpsilon(thisIP)
375 <          ShapeMap%Shapes(current)%sigma = getSigma(thisIP)
374 >          ShapeMap%Shapes(current)%epsilon = getEpsilon(myATID)
375 >          ShapeMap%Shapes(current)%sigma = getSigma(myATID)
376  
377         endif
378  
# Line 378 | Line 380 | contains  
380  
381      haveShapeMap = .true.
382  
383 + !    do i = 1, ShapeMap%n_shapes
384 + !       write(*,*) 'i = ', i, ' isLJ = ', ShapeMap%Shapes(i)%isLJ
385 + !    end do
386 +
387    end subroutine complete_Shape_FF
388  
389 +  function getShapeCut(atomID) result(cutValue)
390 +    integer, intent(in) :: atomID
391 +    real(kind=dp) :: cutValue, whoopdedoo
392 +
393 +    !! this is just a placeholder for a cutoff value, hopefully we'll
394 +    !! develop a method to calculate a sensible value
395 +    whoopdedoo = 9.0_dp
396 +
397 +    cutValue = whoopdedoo
398 +
399 +  end function getShapeCut
400 +
401    subroutine do_shape_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
402         pot, A, f, t, do_pot)
403  
# Line 479 | Line 497 | contains  
497      real (kind=dp) :: fxii, fyii, fzii, fxij, fyij, fzij
498      real (kind=dp) :: fxji, fyji, fzji, fxjj, fyjj, fzjj
499      real (kind=dp) :: fxradial, fyradial, fzradial
500 +
501 +    real (kind=dp) :: xihat, yihat, zihat, xjhat, yjhat, zjhat
502  
503      real (kind=dp) :: plm_i(0:LMAX,0:MMAX), dlm_i(0:LMAX,0:MMAX)
504      real (kind=dp) :: plm_j(0:LMAX,0:MMAX), dlm_j(0:LMAX,0:MMAX)
# Line 492 | Line 512 | contains  
512  
513      !! We assume that the rotation matrices have already been calculated
514      !! and placed in the A array.
495 –
515      r3 = r2*rij
516      r5 = r3*r2
517  
518      drdxi = -d(1) / rij
519      drdyi = -d(2) / rij
520      drdzi = -d(3) / rij
521 +    drduxi = 0.0_dp
522 +    drduyi = 0.0_dp
523 +    drduzi = 0.0_dp
524  
525      drdxj = d(1) / rij
526      drdyj = d(2) / rij
527      drdzj = d(3) / rij
528 +    drduxj = 0.0_dp
529 +    drduyj = 0.0_dp
530 +    drduzj = 0.0_dp
531  
532      ! find the atom type id (atid) for each atom:
533   #ifdef IS_MPI
# Line 516 | Line 541 | contains  
541      ! use the atid to find the shape type (st) for each atom:
542      st1 = ShapeMap%atidToShape(atid1)
543      st2 = ShapeMap%atidToShape(atid2)
544 +    
545 + !    write(*,*) atom1, atom2, atid1, atid2, st1, st2, ShapeMap%Shapes(st1)%isLJ, ShapeMap%Shapes(st2)%isLJ
546  
547      if (ShapeMap%Shapes(st1)%isLJ) then
548  
549         sigma_i = ShapeMap%Shapes(st1)%sigma
550         s_i = ShapeMap%Shapes(st1)%sigma
551         eps_i = ShapeMap%Shapes(st1)%epsilon
552 <       dsigmaidx = 0.0d0
553 <       dsigmaidy = 0.0d0
554 <       dsigmaidz = 0.0d0
555 <       dsigmaidux = 0.0d0
556 <       dsigmaiduy = 0.0d0
557 <       dsigmaiduz = 0.0d0
558 <       dsidx = 0.0d0
559 <       dsidy = 0.0d0
560 <       dsidz = 0.0d0
561 <       dsidux = 0.0d0
562 <       dsiduy = 0.0d0
563 <       dsiduz = 0.0d0
564 <       depsidx = 0.0d0
565 <       depsidy = 0.0d0
566 <       depsidz = 0.0d0
567 <       depsidux = 0.0d0
568 <       depsiduy = 0.0d0
569 <       depsiduz = 0.0d0
552 >       dsigmaidx = 0.0_dp
553 >       dsigmaidy = 0.0_dp
554 >       dsigmaidz = 0.0_dp
555 >       dsigmaidux = 0.0_dp
556 >       dsigmaiduy = 0.0_dp
557 >       dsigmaiduz = 0.0_dp
558 >       dsidx = 0.0_dp
559 >       dsidy = 0.0_dp
560 >       dsidz = 0.0_dp
561 >       dsidux = 0.0_dp
562 >       dsiduy = 0.0_dp
563 >       dsiduz = 0.0_dp
564 >       depsidx = 0.0_dp
565 >       depsidy = 0.0_dp
566 >       depsidz = 0.0_dp
567 >       depsidux = 0.0_dp
568 >       depsiduy = 0.0_dp
569 >       depsiduz = 0.0_dp
570      else
571  
572   #ifdef IS_MPI
# Line 559 | Line 586 | contains  
586         zi = a(7,atom1)*d(1) + a(8,atom1)*d(2) + a(9,atom1)*d(3)
587  
588   #endif
589 <
589 >       xihat = xi / rij
590 >       yihat = yi / rij
591 >       zihat = zi / rij
592         xi2 = xi*xi
593         yi2 = yi*yi
594         zi2 = zi*zi            
# Line 570 | Line 599 | contains  
599  
600         dctidx = - zi * xi / r3
601         dctidy = - zi * yi / r3
602 <       dctidz = 1.0d0 / rij - zi2 / r3
603 <       dctidux = - (zi * xi2) / r3
604 <       dctiduy = - (zi * yi2) / r3
605 <       dctiduz = zi / rij - (zi2 * zi) / r3
602 >       dctidz = 1.0_dp / rij - zi2 / r3
603 >       dctidux = yi / rij ! - (zi * xi2) / r3
604 >       dctiduy = -xi / rij !- (zi * yi2) / r3
605 >       dctiduz = 0.0_dp !zi / rij - (zi2 * zi) / r3
606  
607         ! this is an attempt to try to truncate the singularity when
608         ! sin(theta) is near 0.0:
609  
610         sti2 = 1.0_dp - cti*cti
611 <       if (dabs(sti2) .lt. 1.0d-12) then
612 <          proji = sqrt(rij * 1.0d-12)
613 <          dcpidx = 1.0d0 / proji
614 <          dcpidy = 0.0d0
611 >       if (abs(sti2) .lt. 1.0e-12_dp) then
612 >          proji = sqrt(rij * 1.0e-12_dp)
613 >          dcpidx = 1.0_dp / proji
614 >          dcpidy = 0.0_dp
615            dcpidux = xi / proji
616 <          dcpiduy = 0.0d0
617 <          dspidx = 0.0d0
618 <          dspidy = 1.0d0 / proji
619 <          dspidux = 0.0d0
616 >          dcpiduy = 0.0_dp
617 >          dspidx = 0.0_dp
618 >          dspidy = 1.0_dp / proji
619 >          dspidux = 0.0_dp
620            dspiduy = yi / proji
621         else
622            proji = sqrt(xi2 + yi2)
623            proji3 = proji*proji*proji
624 <          dcpidx = 1.0d0 / proji - xi2 / proji3
624 >          dcpidx = 1.0_dp / proji - xi2 / proji3
625            dcpidy = - xi * yi / proji3
626            dcpidux = xi / proji - (xi2 * xi) / proji3
627            dcpiduy = - (xi * yi2) / proji3
628            dspidx = - xi * yi / proji3
629 <          dspidy = 1.0d0 / proji - yi2 / proji3
629 >          dspidy = 1.0_dp / proji - yi2 / proji3
630            dspidux = - (yi * xi2) / proji3
631            dspiduy = yi / proji - (yi2 * yi) / proji3
632         endif
633  
634         cpi = xi / proji
635 <       dcpidz = 0.0d0
636 <       dcpiduz = 0.0d0
635 >       dcpidz = 0.0_dp
636 >       dcpiduz = 0.0_dp
637  
638         spi = yi / proji
639 <       dspidz = 0.0d0
640 <       dspiduz = 0.0d0
639 >       dspidz = 0.0_dp
640 >       dspiduz = 0.0_dp
641  
642         call Associated_Legendre(cti, ShapeMap%Shapes(st1)%bigM, &
643              ShapeMap%Shapes(st1)%bigL, LMAX, &
# Line 619 | Line 648 | contains  
648         call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, MMAX, &
649              CHEBYSHEV_UN, um_i, dum_i)
650  
651 <       sigma_i = 0.0d0
652 <       s_i = 0.0d0
653 <       eps_i = 0.0d0
654 <       dsigmaidx = 0.0d0
655 <       dsigmaidy = 0.0d0
656 <       dsigmaidz = 0.0d0
657 <       dsigmaidux = 0.0d0
658 <       dsigmaiduy = 0.0d0
659 <       dsigmaiduz = 0.0d0
660 <       dsidx = 0.0d0
661 <       dsidy = 0.0d0
662 <       dsidz = 0.0d0
663 <       dsidux = 0.0d0
664 <       dsiduy = 0.0d0
665 <       dsiduz = 0.0d0
666 <       depsidx = 0.0d0
667 <       depsidy = 0.0d0
668 <       depsidz = 0.0d0
669 <       depsidux = 0.0d0
670 <       depsiduy = 0.0d0
671 <       depsiduz = 0.0d0
651 >       sigma_i = 0.0_dp
652 >       s_i = 0.0_dp
653 >       eps_i = 0.0_dp
654 >       dsigmaidx = 0.0_dp
655 >       dsigmaidy = 0.0_dp
656 >       dsigmaidz = 0.0_dp
657 >       dsigmaidux = 0.0_dp
658 >       dsigmaiduy = 0.0_dp
659 >       dsigmaiduz = 0.0_dp
660 >       dsidx = 0.0_dp
661 >       dsidy = 0.0_dp
662 >       dsidz = 0.0_dp
663 >       dsidux = 0.0_dp
664 >       dsiduy = 0.0_dp
665 >       dsiduz = 0.0_dp
666 >       depsidx = 0.0_dp
667 >       depsidy = 0.0_dp
668 >       depsidz = 0.0_dp
669 >       depsidux = 0.0_dp
670 >       depsiduy = 0.0_dp
671 >       depsiduz = 0.0_dp
672  
673         do lm = 1, ShapeMap%Shapes(st1)%nContactFuncs
674            l = ShapeMap%Shapes(st1)%ContactFuncLValue(lm)
# Line 652 | Line 681 | contains  
681               dPhuncdX = coeff * dtm_i(m) * dcpidx
682               dPhuncdY = coeff * dtm_i(m) * dcpidy
683               dPhuncdZ = coeff * dtm_i(m) * dcpidz
684 <             dPhuncdUz = coeff * dtm_i(m) * dcpidux
684 >             dPhuncdUx = coeff * dtm_i(m) * dcpidux
685               dPhuncdUy = coeff * dtm_i(m) * dcpiduy
686               dPhuncdUz = coeff * dtm_i(m) * dcpiduz
687            else
# Line 666 | Line 695 | contains  
695            endif
696  
697            sigma_i = sigma_i + plm_i(m,l)*Phunc
698 <
698 > !!$          write(*,*) 'dsigmaidux = ', dsigmaidux
699 > !!$          write(*,*) 'Phunc = ', Phunc
700            dsigmaidx = dsigmaidx + plm_i(m,l)*dPhuncdX + &
701                 Phunc * dlm_i(m,l) * dctidx
702            dsigmaidy = dsigmaidy + plm_i(m,l)*dPhuncdY + &
703                 Phunc * dlm_i(m,l) * dctidy
704            dsigmaidz = dsigmaidz + plm_i(m,l)*dPhuncdZ + &
705                 Phunc * dlm_i(m,l) * dctidz
676 –
706            dsigmaidux = dsigmaidux + plm_i(m,l)* dPhuncdUx + &
707                 Phunc * dlm_i(m,l) * dctidux
708            dsigmaiduy = dsigmaiduy + plm_i(m,l)* dPhuncdUy + &
709                 Phunc * dlm_i(m,l) * dctiduy
710            dsigmaiduz = dsigmaiduz + plm_i(m,l)* dPhuncdUz + &
711                 Phunc * dlm_i(m,l) * dctiduz
712 <
712 > !!$          write(*,*) 'dsigmaidux = ', dsigmaidux, '; dPhuncdUx = ', dPhuncdUx, &
713 > !!$                     '; dctidux = ', dctidux, '; plm_i(m,l) = ', plm_i(m,l), &
714 > !!$                     '; dlm_i(m,l) = ', dlm_i(m,l), '; m = ', m, '; l = ', l
715         end do
716  
717         do lm = 1, ShapeMap%Shapes(st1)%nRangeFuncs
# Line 694 | Line 725 | contains  
725               dPhuncdX = coeff * dtm_i(m) * dcpidx
726               dPhuncdY = coeff * dtm_i(m) * dcpidy
727               dPhuncdZ = coeff * dtm_i(m) * dcpidz
728 <             dPhuncdUz = coeff * dtm_i(m) * dcpidux
728 >             dPhuncdUx = coeff * dtm_i(m) * dcpidux
729               dPhuncdUy = coeff * dtm_i(m) * dcpiduy
730               dPhuncdUz = coeff * dtm_i(m) * dcpiduz
731            else
# Line 736 | Line 767 | contains  
767               dPhuncdX = coeff * dtm_i(m) * dcpidx
768               dPhuncdY = coeff * dtm_i(m) * dcpidy
769               dPhuncdZ = coeff * dtm_i(m) * dcpidz
770 <             dPhuncdUz = coeff * dtm_i(m) * dcpidux
770 >             dPhuncdUx = coeff * dtm_i(m) * dcpidux
771               dPhuncdUy = coeff * dtm_i(m) * dcpiduy
772               dPhuncdUz = coeff * dtm_i(m) * dcpiduz
773            else
# Line 775 | Line 806 | contains  
806         sigma_j = ShapeMap%Shapes(st2)%sigma
807         s_j = ShapeMap%Shapes(st2)%sigma
808         eps_j = ShapeMap%Shapes(st2)%epsilon
809 <       dsigmajdx = 0.0d0
810 <       dsigmajdy = 0.0d0
811 <       dsigmajdz = 0.0d0
812 <       dsigmajdux = 0.0d0
813 <       dsigmajduy = 0.0d0
814 <       dsigmajduz = 0.0d0
815 <       dsjdx = 0.0d0
816 <       dsjdy = 0.0d0
817 <       dsjdz = 0.0d0
818 <       dsjdux = 0.0d0
819 <       dsjduy = 0.0d0
820 <       dsjduz = 0.0d0
821 <       depsjdx = 0.0d0
822 <       depsjdy = 0.0d0
823 <       depsjdz = 0.0d0
824 <       depsjdux = 0.0d0
825 <       depsjduy = 0.0d0
826 <       depsjduz = 0.0d0
809 >       dsigmajdx = 0.0_dp
810 >       dsigmajdy = 0.0_dp
811 >       dsigmajdz = 0.0_dp
812 >       dsigmajdux = 0.0_dp
813 >       dsigmajduy = 0.0_dp
814 >       dsigmajduz = 0.0_dp
815 >       dsjdx = 0.0_dp
816 >       dsjdy = 0.0_dp
817 >       dsjdz = 0.0_dp
818 >       dsjdux = 0.0_dp
819 >       dsjduy = 0.0_dp
820 >       dsjduz = 0.0_dp
821 >       depsjdx = 0.0_dp
822 >       depsjdy = 0.0_dp
823 >       depsjdz = 0.0_dp
824 >       depsjdux = 0.0_dp
825 >       depsjduy = 0.0_dp
826 >       depsjduz = 0.0_dp
827      else
828  
829   #ifdef IS_MPI
# Line 813 | Line 844 | contains  
844         zj = -(a(7,atom2)*d(1) + a(8,atom2)*d(2) + a(9,atom2)*d(3))
845   #endif
846  
847 +       xjhat = xj / rij
848 +       yjhat = yj / rij
849 +       zjhat = zj / rij
850         xj2 = xj*xj
851         yj2 = yj*yj
852         zj2 = zj*zj
# Line 823 | Line 857 | contains  
857  
858         dctjdx = - zj * xj / r3
859         dctjdy = - zj * yj / r3
860 <       dctjdz = 1.0d0 / rij - zj2 / r3
861 <       dctjdux = - (zi * xj2) / r3
862 <       dctjduy = - (zj * yj2) / r3
863 <       dctjduz = zj / rij - (zj2 * zj) / r3
860 >       dctjdz = 1.0_dp / rij - zj2 / r3
861 >       dctjdux = yj / rij !- (zi * xj2) / r3
862 >       dctjduy = -xj / rij !- (zj * yj2) / r3
863 >       dctjduz = 0.0_dp !zj / rij - (zj2 * zj) / r3
864  
865         ! this is an attempt to try to truncate the singularity when
866         ! sin(theta) is near 0.0:
867  
868         stj2 = 1.0_dp - ctj*ctj
869 <       if (dabs(stj2) .lt. 1.0d-12) then
870 <          projj = sqrt(rij * 1.0d-12)
871 <          dcpjdx = 1.0d0 / projj
872 <          dcpjdy = 0.0d0
869 >       if (abs(stj2) .lt. 1.0e-12_dp) then
870 >          projj = sqrt(rij * 1.0e-12_dp)
871 >          dcpjdx = 1.0_dp / projj
872 >          dcpjdy = 0.0_dp
873            dcpjdux = xj / projj
874 <          dcpjduy = 0.0d0
875 <          dspjdx = 0.0d0
876 <          dspjdy = 1.0d0 / projj
877 <          dspjdux = 0.0d0
874 >          dcpjduy = 0.0_dp
875 >          dspjdx = 0.0_dp
876 >          dspjdy = 1.0_dp / projj
877 >          dspjdux = 0.0_dp
878            dspjduy = yj / projj
879         else
880            projj = sqrt(xj2 + yj2)
881            projj3 = projj*projj*projj
882 <          dcpjdx = 1.0d0 / projj - xj2 / projj3
882 >          dcpjdx = 1.0_dp / projj - xj2 / projj3
883            dcpjdy = - xj * yj / projj3
884            dcpjdux = xj / projj - (xj2 * xj) / projj3
885            dcpjduy = - (xj * yj2) / projj3
886            dspjdx = - xj * yj / projj3
887 <          dspjdy = 1.0d0 / projj - yj2 / projj3
887 >          dspjdy = 1.0_dp / projj - yj2 / projj3
888            dspjdux = - (yj * xj2) / projj3
889            dspjduy = yj / projj - (yj2 * yj) / projj3
890         endif
891  
892         cpj = xj / projj
893 <       dcpjdz = 0.0d0
894 <       dcpjduz = 0.0d0
893 >       dcpjdz = 0.0_dp
894 >       dcpjduz = 0.0_dp
895  
896         spj = yj / projj
897 <       dspjdz = 0.0d0
898 <       dspjduz = 0.0d0
897 >       dspjdz = 0.0_dp
898 >       dspjduz = 0.0_dp
899  
900  
901 <       write(*,*) 'dcpdu = ' ,dcpidux, dcpiduy, dcpiduz
902 <       write(*,*) 'dcpdu = ' ,dcpjdux, dcpjduy, dcpjduz
901 > !       write(*,*) 'dcpdu = ' ,dcpidux, dcpiduy, dcpiduz
902 > !       write(*,*) 'dcpdu = ' ,dcpjdux, dcpjduy, dcpjduz
903         call Associated_Legendre(ctj, ShapeMap%Shapes(st2)%bigM, &
904              ShapeMap%Shapes(st2)%bigL, LMAX, &
905              plm_j, dlm_j)
# Line 875 | Line 909 | contains  
909         call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, MMAX, &
910              CHEBYSHEV_UN, um_j, dum_j)
911  
912 <       sigma_j = 0.0d0
913 <       s_j = 0.0d0
914 <       eps_j = 0.0d0
915 <       dsigmajdx = 0.0d0
916 <       dsigmajdy = 0.0d0
917 <       dsigmajdz = 0.0d0
918 <       dsigmajdux = 0.0d0
919 <       dsigmajduy = 0.0d0
920 <       dsigmajduz = 0.0d0
921 <       dsjdx = 0.0d0
922 <       dsjdy = 0.0d0
923 <       dsjdz = 0.0d0
924 <       dsjdux = 0.0d0
925 <       dsjduy = 0.0d0
926 <       dsjduz = 0.0d0
927 <       depsjdx = 0.0d0
928 <       depsjdy = 0.0d0
929 <       depsjdz = 0.0d0
930 <       depsjdux = 0.0d0
931 <       depsjduy = 0.0d0
932 <       depsjduz = 0.0d0
912 >       sigma_j = 0.0_dp
913 >       s_j = 0.0_dp
914 >       eps_j = 0.0_dp
915 >       dsigmajdx = 0.0_dp
916 >       dsigmajdy = 0.0_dp
917 >       dsigmajdz = 0.0_dp
918 >       dsigmajdux = 0.0_dp
919 >       dsigmajduy = 0.0_dp
920 >       dsigmajduz = 0.0_dp
921 >       dsjdx = 0.0_dp
922 >       dsjdy = 0.0_dp
923 >       dsjdz = 0.0_dp
924 >       dsjdux = 0.0_dp
925 >       dsjduy = 0.0_dp
926 >       dsjduz = 0.0_dp
927 >       depsjdx = 0.0_dp
928 >       depsjdy = 0.0_dp
929 >       depsjdz = 0.0_dp
930 >       depsjdux = 0.0_dp
931 >       depsjduy = 0.0_dp
932 >       depsjduz = 0.0_dp
933  
934         do lm = 1, ShapeMap%Shapes(st2)%nContactFuncs
935            l = ShapeMap%Shapes(st2)%ContactFuncLValue(lm)
# Line 908 | Line 942 | contains  
942               dPhuncdX = coeff * dtm_j(m) * dcpjdx
943               dPhuncdY = coeff * dtm_j(m) * dcpjdy
944               dPhuncdZ = coeff * dtm_j(m) * dcpjdz
945 <             dPhuncdUz = coeff * dtm_j(m) * dcpjdux
945 >             dPhuncdUx = coeff * dtm_j(m) * dcpjdux
946               dPhuncdUy = coeff * dtm_j(m) * dcpjduy
947               dPhuncdUz = coeff * dtm_j(m) * dcpjduz
948            else
# Line 950 | Line 984 | contains  
984               dPhuncdX = coeff * dtm_j(m) * dcpjdx
985               dPhuncdY = coeff * dtm_j(m) * dcpjdy
986               dPhuncdZ = coeff * dtm_j(m) * dcpjdz
987 <             dPhuncdUz = coeff * dtm_j(m) * dcpjdux
987 >             dPhuncdUx = coeff * dtm_j(m) * dcpjdux
988               dPhuncdUy = coeff * dtm_j(m) * dcpjduy
989               dPhuncdUz = coeff * dtm_j(m) * dcpjduz
990            else
# Line 1005 | Line 1039 | contains  
1039               dPhuncdUz = coeff*(spj * dum_j(m-1)*dcpjduz + dspjduz *um_j(m-1))
1040            endif
1041  
1042 <          write(*,*) 'l,m = ', l, m, coeff, dPhuncdUx, dPhuncdUy, dPhuncdUz
1042 > !          write(*,*) 'l,m = ', l, m, coeff, dPhuncdUx, dPhuncdUy, dPhuncdUz
1043  
1044            eps_j = eps_j + plm_j(m,l)*Phunc
1045  
# Line 1030 | Line 1064 | contains  
1064      ! phew, now let's assemble the potential energy:
1065  
1066      sigma = 0.5*(sigma_i + sigma_j)
1067 <
1067 > !    write(*,*) sigma_i, ' = sigma_i; ', sigma_j, ' = sigma_j'
1068      dsigmadxi = 0.5*dsigmaidx
1069      dsigmadyi = 0.5*dsigmaidy
1070      dsigmadzi = 0.5*dsigmaidz
# Line 1062 | Line 1096 | contains  
1096      dsduzj = 0.5*dsjduz
1097  
1098      eps = sqrt(eps_i * eps_j)
1099 <
1100 <    depsdxi = eps_j * depsidx / (2.0d0 * eps)
1101 <    depsdyi = eps_j * depsidy / (2.0d0 * eps)
1102 <    depsdzi = eps_j * depsidz / (2.0d0 * eps)
1103 <    depsduxi = eps_j * depsidux / (2.0d0 * eps)
1104 <    depsduyi = eps_j * depsiduy / (2.0d0 * eps)
1105 <    depsduzi = eps_j * depsiduz / (2.0d0 * eps)
1099 > !!$    write(*,*) 'dsidu = ', dsidux, dsiduy, dsiduz
1100 > !!$    write(*,*) 'dsigidu = ', dsigmaidux, dsigmaiduy, dsigmaiduz
1101 > !!$    write(*,*) sigma_j, ' is sigma j; ', s_j, ' is s j; ', eps_j, ' is eps j'
1102 >    depsdxi = eps_j * depsidx / (2.0_dp * eps)
1103 >    depsdyi = eps_j * depsidy / (2.0_dp * eps)
1104 >    depsdzi = eps_j * depsidz / (2.0_dp * eps)
1105 >    depsduxi = eps_j * depsidux / (2.0_dp * eps)
1106 >    depsduyi = eps_j * depsiduy / (2.0_dp * eps)
1107 >    depsduzi = eps_j * depsiduz / (2.0_dp * eps)
1108  
1109 <    depsdxj = eps_i * depsjdx / (2.0d0 * eps)
1110 <    depsdyj = eps_i * depsjdy / (2.0d0 * eps)
1111 <    depsdzj = eps_i * depsjdz / (2.0d0 * eps)
1112 <    depsduxj = eps_i * depsjdux / (2.0d0 * eps)
1113 <    depsduyj = eps_i * depsjduy / (2.0d0 * eps)
1114 <    depsduzj = eps_i * depsjduz / (2.0d0 * eps)
1109 >    depsdxj = eps_i * depsjdx / (2.0_dp * eps)
1110 >    depsdyj = eps_i * depsjdy / (2.0_dp * eps)
1111 >    depsdzj = eps_i * depsjdz / (2.0_dp * eps)
1112 >    depsduxj = eps_i * depsjdux / (2.0_dp * eps)
1113 >    depsduyj = eps_i * depsjduy / (2.0_dp * eps)
1114 >    depsduzj = eps_i * depsjduz / (2.0_dp * eps)
1115  
1116   !!$    write(*,*) 'depsidu = ', depsidux, depsiduy, depsiduz
1117 +
1118   !!$    write(*,*) 'depsjdu = ', depsjdux, depsjduy, depsjduz
1082 – !!$
1083 – !!$    write(*,*) 'depsdui = ', depsduxi, depsduyi, depsduzi
1119   !!$    write(*,*) 'depsduj = ', depsduxj, depsduyj, depsduzj
1120   !!$
1121   !!$    write(*,*) 's, sig, eps = ', s, sigma, eps
# Line 1088 | Line 1123 | contains  
1123      rtdenom = rij-sigma+s
1124      rt = s / rtdenom
1125  
1126 <    drtdxi = (dsdxi + rt * (drdxi - dsigmadxi + dsdxi)) / rtdenom
1127 <    drtdyi = (dsdyi + rt * (drdyi - dsigmadyi + dsdyi)) / rtdenom
1128 <    drtdzi = (dsdzi + rt * (drdzi - dsigmadzi + dsdzi)) / rtdenom
1129 <    drtduxi = (dsduxi + rt * (drduxi - dsigmaduxi + dsduxi)) / rtdenom
1130 <    drtduyi = (dsduyi + rt * (drduyi - dsigmaduyi + dsduyi)) / rtdenom
1131 <    drtduzi = (dsduzi + rt * (drduzi - dsigmaduzi + dsduzi)) / rtdenom
1132 <    drtdxj = (dsdxj + rt * (drdxj - dsigmadxj + dsdxj)) / rtdenom
1133 <    drtdyj = (dsdyj + rt * (drdyj - dsigmadyj + dsdyj)) / rtdenom
1134 <    drtdzj = (dsdzj + rt * (drdzj - dsigmadzj + dsdzj)) / rtdenom
1135 <    drtduxj = (dsduxj + rt * (drduxj - dsigmaduxj + dsduxj)) / rtdenom
1136 <    drtduyj = (dsduyj + rt * (drduyj - dsigmaduyj + dsduyj)) / rtdenom
1137 <    drtduzj = (dsduzj + rt * (drduzj - dsigmaduzj + dsduzj)) / rtdenom
1126 >    drtdxi = (dsdxi - rt * (drdxi - dsigmadxi + dsdxi)) / rtdenom
1127 >    drtdyi = (dsdyi - rt * (drdyi - dsigmadyi + dsdyi)) / rtdenom
1128 >    drtdzi = (dsdzi - rt * (drdzi - dsigmadzi + dsdzi)) / rtdenom
1129 >    drtduxi = (dsduxi - rt * (drduxi - dsigmaduxi + dsduxi)) / rtdenom
1130 >    drtduyi = (dsduyi - rt * (drduyi - dsigmaduyi + dsduyi)) / rtdenom
1131 >    drtduzi = (dsduzi - rt * (drduzi - dsigmaduzi + dsduzi)) / rtdenom
1132 >    drtdxj = (dsdxj - rt * (drdxj - dsigmadxj + dsdxj)) / rtdenom
1133 >    drtdyj = (dsdyj - rt * (drdyj - dsigmadyj + dsdyj)) / rtdenom
1134 >    drtdzj = (dsdzj - rt * (drdzj - dsigmadzj + dsdzj)) / rtdenom
1135 >    drtduxj = (dsduxj - rt * (drduxj - dsigmaduxj + dsduxj)) / rtdenom
1136 >    drtduyj = (dsduyj - rt * (drduyj - dsigmaduyj + dsduyj)) / rtdenom
1137 >    drtduzj = (dsduzj - rt * (drduzj - dsigmaduzj + dsduzj)) / rtdenom
1138  
1139 + !!$    write(*,*) 'drtd_i = ', drtdxi, drtdyi, drtdzi
1140 + !!$    write(*,*) 'drtdu_j = ', drtduxj, drtduyj, drtduzj
1141 +
1142      rt2 = rt*rt
1143      rt3 = rt2*rt
1144      rt5 = rt2*rt3
# Line 1109 | Line 1147 | contains  
1147      rt12 = rt6*rt6
1148      rt126 = rt12 - rt6
1149  
1150 <    pot_temp = 4.0d0 * eps * rt126
1150 >    pot_temp = 4.0_dp * eps * rt126
1151  
1152      vpair = vpair + pot_temp
1153      if (do_pot) then
1154   #ifdef IS_MPI
1155 <       pot_row(atom1) = pot_row(atom1) + 0.5d0*pot_temp*sw
1156 <       pot_col(atom2) = pot_col(atom2) + 0.5d0*pot_temp*sw
1155 >       pot_row(VDW_POT,atom1) = pot_row(VDW_POT,atom1) + 0.5_dp*pot_temp*sw
1156 >       pot_col(VDW_POT,atom2) = pot_col(VDW_POT,atom2) + 0.5_dp*pot_temp*sw
1157   #else
1158         pot = pot + pot_temp*sw
1159   #endif
# Line 1123 | Line 1161 | contains  
1161  
1162   !!$    write(*,*) 'drtdu, depsdu = ', drtduxi, depsduxi
1163  
1164 <    dvdxi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdxi + 4.0d0*depsdxi*rt126
1165 <    dvdyi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdyi + 4.0d0*depsdyi*rt126
1166 <    dvdzi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdzi + 4.0d0*depsdzi*rt126
1167 <    dvduxi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduxi + 4.0d0*depsduxi*rt126
1168 <    dvduyi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduyi + 4.0d0*depsduyi*rt126
1169 <    dvduzi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduzi + 4.0d0*depsduzi*rt126
1132 <
1133 <    dvdxj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdxj + 4.0d0*depsdxj*rt126
1134 <    dvdyj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdyj + 4.0d0*depsdyj*rt126
1135 <    dvdzj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdzj + 4.0d0*depsdzj*rt126
1136 <    dvduxj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduxj + 4.0d0*depsduxj*rt126
1137 <    dvduyj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduyj + 4.0d0*depsduyj*rt126
1138 <    dvduzj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduzj + 4.0d0*depsduzj*rt126
1164 >    dvdxi = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtdxi + 4.0_dp*depsdxi*rt126
1165 >    dvdyi = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtdyi + 4.0_dp*depsdyi*rt126
1166 >    dvdzi = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtdzi + 4.0_dp*depsdzi*rt126
1167 >    dvduxi = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtduxi + 4.0_dp*depsduxi*rt126
1168 >    dvduyi = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtduyi + 4.0_dp*depsduyi*rt126
1169 >    dvduzi = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtduzi + 4.0_dp*depsduzi*rt126
1170  
1171 +    dvdxj = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtdxj + 4.0_dp*depsdxj*rt126
1172 +    dvdyj = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtdyj + 4.0_dp*depsdyj*rt126
1173 +    dvdzj = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtdzj + 4.0_dp*depsdzj*rt126
1174 +    dvduxj = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtduxj + 4.0_dp*depsduxj*rt126
1175 +    dvduyj = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtduyj + 4.0_dp*depsduyj*rt126
1176 +    dvduzj = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtduzj + 4.0_dp*depsduzj*rt126
1177 + !!$    write(*,*) 'drtduxi = ', drtduxi, ' depsduxi = ', depsduxi
1178      ! do the torques first since they are easy:
1179      ! remember that these are still in the body fixed axes
1180  
1181 +    txi = 0.0_dp
1182 +    tyi = 0.0_dp
1183 +    tzi = 0.0_dp
1184  
1185 < !!$    write(*,*) 'sw = ', sw
1186 < !!$    write(*,*) 'dvdu1 = ', dvduxi, dvduyi, dvduzi
1187 < !!$    write(*,*) 'dvdu2 = ', dvduxj, dvduyj, dvduzj
1147 < !!$
1148 <    txi =  (dvduzi - dvduyi) * sw
1149 <    tyi =  (dvduxi - dvduzi) * sw
1150 <    tzi =  (dvduyi - dvduxi) * sw
1185 >    txj = 0.0_dp
1186 >    tyj = 0.0_dp
1187 >    tzj = 0.0_dp
1188  
1189 <    txj = (dvduzj - dvduyj) * sw
1190 <    tyj = (dvduxj - dvduzj) * sw
1191 <    tzj = (dvduyj - dvduxj) * sw
1189 >    txi = (dvduyi - dvduzi) * sw
1190 >    tyi = (dvduzi - dvduxi) * sw
1191 >    tzi = (dvduxi - dvduyi) * sw
1192  
1193 < !!$    txi = -dvduxi * sw
1194 < !!$    tyi = -dvduyi * sw
1195 < !!$    tzi = -dvduzi * sw
1193 >    txj = (dvduyj - dvduzj) * sw
1194 >    tyj = (dvduzj - dvduxj) * sw
1195 >    tzj = (dvduxj - dvduyj) * sw
1196 >
1197 > !!$    txi = dvduxi * sw
1198 > !!$    tyi = dvduyi * sw
1199 > !!$    tzi = dvduzi * sw
1200   !!$
1201   !!$    txj = dvduxj * sw
1202   !!$    tyj = dvduyj * sw
# Line 1188 | Line 1229 | contains  
1229      t(1,atom2) = t(1,atom2) + a(1,atom2)*txj + a(4,atom2)*tyj + a(7,atom2)*tzj
1230      t(2,atom2) = t(2,atom2) + a(2,atom2)*txj + a(5,atom2)*tyj + a(8,atom2)*tzj
1231      t(3,atom2) = t(3,atom2) + a(3,atom2)*txj + a(6,atom2)*tyj + a(9,atom2)*tzj
1232 +
1233   #endif
1234      ! Now, on to the forces:
1235  
1236      ! first rotate the i terms back into the lab frame:
1237  
1238 <    fxi = dvdxi * sw
1239 <    fyi = dvdyi * sw
1240 <    fzi = dvdzi * sw
1238 >    fxi = -dvdxi * sw
1239 >    fyi = -dvdyi * sw
1240 >    fzi = -dvdzi * sw
1241  
1242 <    fxj = dvdxj * sw
1243 <    fyj = dvdyj * sw
1244 <    fzj = dvdzj * sw
1242 >    fxj = -dvdxj * sw
1243 >    fyj = -dvdyj * sw
1244 >    fzj = -dvdzj * sw
1245  
1246 +
1247   #ifdef IS_MPI
1248      fxii = a_Row(1,atom1)*fxi + a_Row(4,atom1)*fyi + a_Row(7,atom1)*fzi
1249      fyii = a_Row(2,atom1)*fxi + a_Row(5,atom1)*fyi + a_Row(8,atom1)*fzi
# Line 1227 | Line 1270 | contains  
1270      fyji = -fyjj
1271      fzji = -fzjj
1272  
1273 <    fxradial = 0.5_dp * (fxii + fxji)
1274 <    fyradial = 0.5_dp * (fyii + fyji)
1275 <    fzradial = 0.5_dp * (fzii + fzji)
1276 <
1273 >    fxradial = (fxii + fxji)
1274 >    fyradial = (fyii + fyji)
1275 >    fzradial = (fzii + fzji)
1276 > !!$    write(*,*) fxradial, ' is fxrad; ', fyradial, ' is fyrad; ', fzradial, 'is fzrad'
1277   #ifdef IS_MPI
1278      f_Row(1,atom1) = f_Row(1,atom1) + fxradial
1279      f_Row(2,atom1) = f_Row(2,atom1) + fyradial
# Line 1300 | Line 1343 | contains  
1343      end DO
1344  
1345      ! start with 0,0:
1346 <    PLM(0,0) = 1.0D0
1346 >    PLM(0,0) = 1.0_DP
1347  
1348      ! x = +/- 1 functions are easy:
1349 <    IF (abs(X).EQ.1.0D0) THEN
1349 >    IF (abs(X).EQ.1.0_DP) THEN
1350         DO I = 1, m
1351            PLM(0, I) = X**I
1352 <          DLM(0, I) = 0.5D0*I*(I+1.0D0)*X**(I+1)
1352 >          DLM(0, I) = 0.5_DP*I*(I+1.0_DP)*X**(I+1)
1353         end DO
1354         DO J = 1, m
1355            DO I = 1, l
1356               IF (I.EQ.1) THEN
1357                  DLM(I, J) = 1.0D+300
1358               ELSE IF (I.EQ.2) THEN
1359 <                DLM(I, J) = -0.25D0*(J+2)*(J+1)*J*(J-1)*X**(J+1)
1359 >                DLM(I, J) = -0.25_DP*(J+2)*(J+1)*J*(J-1)*X**(J+1)
1360               ENDIF
1361            end DO
1362         end DO
# Line 1321 | Line 1364 | contains  
1364      ENDIF
1365  
1366      LS = 1
1367 <    IF (abs(X).GT.1.0D0) LS = -1
1368 <    XQ = sqrt(LS*(1.0D0-X*X))
1369 <    XS = LS*(1.0D0-X*X)
1367 >    IF (abs(X).GT.1.0_DP) LS = -1
1368 >    XQ = sqrt(LS*(1.0_DP-X*X))
1369 >    XS = LS*(1.0_DP-X*X)
1370  
1371      DO I = 1, l
1372 <       PLM(I, I) = -LS*(2.0D0*I-1.0D0)*XQ*PLM(I-1, I-1)
1372 >       PLM(I, I) = -LS*(2.0_DP*I-1.0_DP)*XQ*PLM(I-1, I-1)
1373      enddo
1374  
1375      DO I = 0, l
1376 <       PLM(I, I+1)=(2.0D0*I+1.0D0)*X*PLM(I, I)
1376 >       PLM(I, I+1)=(2.0_DP*I+1.0_DP)*X*PLM(I, I)
1377      enddo
1378  
1379      DO I = 0, l
1380         DO J = I+2, m
1381 <          PLM(I, J)=((2.0D0*J-1.0D0)*X*PLM(I,J-1) - &
1382 <               (I+J-1.0D0)*PLM(I,J-2))/(J-I)
1381 >          PLM(I, J)=((2.0_DP*J-1.0_DP)*X*PLM(I,J-1) - &
1382 >               (I+J-1.0_DP)*PLM(I,J-2))/(J-I)
1383         end DO
1384      end DO
1385  
1386 <    DLM(0, 0)=0.0D0
1386 >    DLM(0, 0)=0.0_DP
1387      DO J = 1, m
1388         DLM(0, J)=LS*J*(PLM(0,J-1)-X*PLM(0,J))/XS
1389      end DO
1390  
1391      DO I = 1, l
1392         DO J = I, m
1393 <          DLM(I,J) = LS*I*X*PLM(I, J)/XS + (J+I)*(J-I+1.0D0)/XQ*PLM(I-1, J)
1393 >          DLM(I,J) = LS*I*X*PLM(I, J)/XS + (J+I)*(J-I+1.0_DP)/XQ*PLM(I-1, J)
1394         end DO
1395      end DO
1396  
# Line 1384 | Line 1427 | contains  
1427      real(kind=8) :: a, b, c, y0, y1, dy0, dy1, yn, dyn
1428      integer :: k
1429  
1430 <    A = 2.0D0
1431 <    B = 0.0D0
1432 <    C = 1.0D0
1433 <    Y0 = 1.0D0
1434 <    Y1 = 2.0D0*X
1435 <    DY0 = 0.0D0
1436 <    DY1 = 2.0D0
1437 <    PL(0) = 1.0D0
1438 <    PL(1) = 2.0D0*X
1439 <    DPL(0) = 0.0D0
1440 <    DPL(1) = 2.0D0
1430 >    A = 2.0_DP
1431 >    B = 0.0_DP
1432 >    C = 1.0_DP
1433 >    Y0 = 1.0_DP
1434 >    Y1 = 2.0_DP*X
1435 >    DY0 = 0.0_DP
1436 >    DY1 = 2.0_DP
1437 >    PL(0) = 1.0_DP
1438 >    PL(1) = 2.0_DP*X
1439 >    DPL(0) = 0.0_DP
1440 >    DPL(1) = 2.0_DP
1441      IF (function_type.EQ.CHEBYSHEV_TN) THEN
1442         Y1 = X
1443 <       DY1 = 1.0D0
1443 >       DY1 = 1.0_DP
1444         PL(1) = X
1445 <       DPL(1) = 1.0D0
1445 >       DPL(1) = 1.0_DP
1446      ELSE IF (function_type.EQ.LAGUERRE) THEN
1447 <       Y1 = 1.0D0-X
1448 <       DY1 = -1.0D0
1449 <       PL(1) = 1.0D0-X
1450 <       DPL(1) = -1.0D0
1447 >       Y1 = 1.0_DP-X
1448 >       DY1 = -1.0_DP
1449 >       PL(1) = 1.0_DP-X
1450 >       DPL(1) = -1.0_DP
1451      ENDIF
1452      DO K = 2, m
1453         IF (function_type.EQ.LAGUERRE) THEN
1454 <          A = -1.0D0/K
1455 <          B = 2.0D0+A
1456 <          C = 1.0D0+A
1454 >          A = -1.0_DP/K
1455 >          B = 2.0_DP+A
1456 >          C = 1.0_DP+A
1457         ELSE IF (function_type.EQ.HERMITE) THEN
1458 <          C = 2.0D0*(K-1.0D0)
1458 >          C = 2.0_DP*(K-1.0_DP)
1459         ENDIF
1460         YN = (A*X+B)*Y1-C*Y0
1461         DYN = A*Y1+(A*X+B)*DY1-C*DY0

Diff Legend

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