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

Comparing trunk/OOPSE-4/src/UseTheForce/DarkSide/electrostatic.F90 (file contents):
Revision 2301 by gezelter, Thu Sep 15 22:05:21 2005 UTC vs.
Revision 2330 by chuckv, Mon Sep 26 17:07:54 2005 UTC

# Line 73 | Line 73 | module electrostatic_module
73  
74    !! variables to handle different summation methods for long-range electrostatics:
75    integer, save :: summationMethod = NONE
76 +  logical, save :: summationMethodChecked = .false.
77    real(kind=DP), save :: defaultCutoff = 0.0_DP
78    logical, save :: haveDefaultCutoff = .false.
79    real(kind=DP), save :: dampingAlpha = 0.0_DP
# Line 82 | Line 83 | module electrostatic_module
83    real(kind=DP), save :: constERFC = 0.0_DP
84    real(kind=DP), save :: constEXP = 0.0_DP
85    logical, save :: haveDWAconstants = .false.
86 +  real(kind=dp), save :: rcuti = 0.0_dp
87 +  real(kind=dp), save :: rcuti2 = 0.0_dp
88 +  real(kind=dp), save :: rcuti3 = 0.0_dp
89 +  real(kind=dp), save :: rcuti4 = 0.0_dp
90 +  logical, save :: is_Undamped_Wolf = .false.
91 +  logical, save :: is_Damped_Wolf = .false.
92 +
93  
94 + ! error function
95 +  double precision, external :: derfc
96  
97    public :: setElectrostaticSummationMethod
98    public :: setElectrostaticCutoffRadius
# Line 117 | Line 127 | contains
127   contains
128  
129    subroutine setElectrostaticSummationMethod(the_ESM)
120 –
130      integer, intent(in) :: the_ESM    
131  
132      if ((the_ESM .le. 0) .or. (the_ESM .gt. REACTION_FIELD)) then
133         call handleError("setElectrostaticSummationMethod", "Unsupported Summation Method")
134      endif
135 +
136 +    summationMethod = the_ESM
137  
138 +    if (summationMethod .eq. UNDAMPED_WOLF) is_Undamped_Wolf = .true.
139 +    if (summationMethod .eq. DAMPED_WOLF) is_Damped_Wolf = .true.
140    end subroutine setElectrostaticSummationMethod
141  
142    subroutine setElectrostaticCutoffRadius(thisRcut)
# Line 355 | Line 368 | contains
368    end function getDipoleMoment
369  
370    subroutine checkSummationMethod()
371 +
372 +    if (.not.haveDefaultCutoff) then
373 +       call handleError("checkSummationMethod", "no Default Cutoff set!")
374 +    endif
375  
376 +    rcuti = 1.0d0 / defaultCutoff
377 +    rcuti2 = rcuti*rcuti
378 +    rcuti3 = rcuti2*rcuti
379 +    rcuti4 = rcuti2*rcuti2
380 +
381      if (summationMethod .eq. DAMPED_WOLF) then
382         if (.not.haveDWAconstants) then
383            
# Line 363 | Line 385 | contains
385               call handleError("checkSummationMethod", "no Damping Alpha set!")
386            endif
387            
388 <          if (.not.have....)
389 <          constEXP =
390 <          constERFC =
388 >          if (.not.haveDefaultCutoff) then
389 >             call handleError("checkSummationMethod", "no Default Cutoff set!")
390 >          endif
391 >
392 >          constEXP = exp(-dampingAlpha*dampingAlpha*defaultCutoff*defaultCutoff)
393 >          constERFC = derfc(dampingAlpha*defaultCutoff)
394            
395            haveDWAconstants = .true.
396         endif
397      endif
398  
399 +    if (summationMethod .eq. REACTION_FIELD) then
400 +       if (.not.haveDielectric) then
401 +          call handleError("checkSummationMethod", "no reaction field Dielectric set!")
402 +       endif
403 +    endif
404 +
405 +    summationMethodChecked = .true.
406    end subroutine checkSummationMethod
407  
408  
409  
410    subroutine doElectrostaticPair(atom1, atom2, d, rij, r2, sw, &
411 <       vpair, fpair, pot, eFrame, f, t, do_pot, corrMethod, rcuti)
411 >       vpair, fpair, pot, eFrame, f, t, do_pot)
412  
413      logical, intent(in) :: do_pot
414  
415      integer, intent(in) :: atom1, atom2
416      integer :: localError
385 –    integer, intent(in) :: corrMethod
417  
418 <    real(kind=dp), intent(in) :: rij, r2, sw, rcuti
418 >    real(kind=dp), intent(in) :: rij, r2, sw
419      real(kind=dp), intent(in), dimension(3) :: d
420      real(kind=dp), intent(inout) :: vpair
421      real(kind=dp), intent(inout), dimension(3) :: fpair
422  
423 <    real( kind = dp ) :: pot, swi
423 >    real( kind = dp ) :: pot
424      real( kind = dp ), dimension(9,nLocal) :: eFrame
425      real( kind = dp ), dimension(3,nLocal) :: f
426      real( kind = dp ), dimension(3,nLocal) :: t
# Line 414 | Line 445 | contains
445      real (kind=dp) :: pref, vterm, epot, dudr, vterm1, vterm2
446      real (kind=dp) :: xhat, yhat, zhat
447      real (kind=dp) :: dudx, dudy, dudz
448 <    real (kind=dp) :: scale, sc2, bigR, switcher, dswitcher
418 <    real (kind=dp) :: rcuti2, rcuti3, rcuti4
448 >    real (kind=dp) :: scale, sc2, bigR
449  
450      if (.not.allocated(ElectrostaticMap)) then
451         call handleError("electrostatic", "no ElectrostaticMap was present before first call of do_electrostatic_pair!")
# Line 424 | Line 454 | contains
454  
455      if (.not.summationMethodChecked) then
456         call checkSummationMethod()
457 +      
458      endif
459  
460  
# Line 443 | Line 474 | contains
474      yhat = d(2) * riji
475      zhat = d(3) * riji
476  
446 –    rcuti2 = rcuti*rcuti
447 –    rcuti3 = rcuti2*rcuti
448 –    rcuti4 = rcuti2*rcuti2
449 –
450 –    swi = 1.0d0 / sw
451 –
477      !! logicals
478      i_is_Charge = ElectrostaticMap(me1)%is_Charge
479      i_is_Dipole = ElectrostaticMap(me1)%is_Dipole
# Line 567 | Line 592 | contains
592         cz_j = uz_j(1)*xhat + uz_j(2)*yhat + uz_j(3)*zhat
593      endif
594    
570 – !!$    switcher = 1.0d0
571 – !!$    dswitcher = 0.0d0
572 – !!$    ebalance = 0.0d0
573 – !!$    ! weaken the dipole interaction at close range for TAP water
574 – !!$    if (j_is_Tap .and. i_is_Tap) then
575 – !!$      call calc_switch(rij, mu_i, switcher, dswitcher)
576 – !!$    endif
577 –
595      epot = 0.0_dp
596      dudx = 0.0_dp
597      dudy = 0.0_dp
# Line 592 | Line 609 | contains
609  
610         if (j_is_Charge) then
611  
612 <          if (corrMethod .eq. 1) then
596 <             vterm = pre11 * q_i * q_j * (riji - rcuti)
612 >          if (summationMethod .eq. UNDAMPED_WOLF) then
613  
614 +             vterm = pre11 * q_i * q_j * (riji - rcuti)
615               vpair = vpair + vterm
616 <             epot = epot + sw * vterm
616 >             epot = epot + sw*vterm
617              
618               dudr  = - sw * pre11 * q_i * q_j * (riji*riji*riji - rcuti2*rcuti)
619              
# Line 605 | Line 622 | contains
622               dudz = dudz + dudr * d(3)
623  
624            else
608 –             vterm = pre11 * q_i * q_j * riji
625  
626 +             vterm = pre11 * q_i * q_j * riji
627               vpair = vpair + vterm
628 <             epot = epot + sw * vterm
628 >             epot = epot + sw*vterm
629              
630               dudr  = - sw * vterm * riji
631              
# Line 622 | Line 639 | contains
639  
640         if (j_is_Dipole) then
641  
642 <          pref = sw * pre12 * q_i * mu_j
642 >          pref = pre12 * q_i * mu_j
643  
644 <          if (corrMethod .eq. 1) then
644 >          if (summationMethod .eq. UNDAMPED_WOLF) then
645               ri2 = riji * riji
646               ri3 = ri2 * riji
647  
648 +             pref = pre12 * q_i * mu_j
649               vterm = - pref * ct_j * (ri2 - rcuti2)
650 <             vpair = vpair + swi*vterm
651 <             epot = epot + vterm
650 >             vpair = vpair + vterm
651 >             epot = epot + sw*vterm
652              
653               !! this has a + sign in the () because the rij vector is
654               !! r_j - r_i and the charge-dipole potential takes the origin
655               !! as the point dipole, which is atom j in this case.
656              
657 <             dudx = dudx - pref * ( ri3*( uz_j(1) - 3.0d0*ct_j*xhat) &
657 >             dudx = dudx - sw*pref * ( ri3*( uz_j(1) - 3.0d0*ct_j*xhat) &
658                    - rcuti3*( uz_j(1) - 3.0d0*ct_j*d(1)*rcuti ) )
659 <             dudy = dudy - pref * ( ri3*( uz_j(2) - 3.0d0*ct_j*yhat) &
659 >             dudy = dudy - sw*pref * ( ri3*( uz_j(2) - 3.0d0*ct_j*yhat) &
660                    - rcuti3*( uz_j(2) - 3.0d0*ct_j*d(2)*rcuti ) )
661 <             dudz = dudz - pref * ( ri3*( uz_j(3) - 3.0d0*ct_j*zhat) &
661 >             dudz = dudz - sw*pref * ( ri3*( uz_j(3) - 3.0d0*ct_j*zhat) &
662                    - rcuti3*( uz_j(3) - 3.0d0*ct_j*d(3)*rcuti ) )
663              
664 <             duduz_j(1) = duduz_j(1) - pref*( ri2*xhat - d(1)*rcuti3 )
665 <             duduz_j(2) = duduz_j(2) - pref*( ri2*yhat - d(2)*rcuti3 )
666 <             duduz_j(3) = duduz_j(3) - pref*( ri2*zhat - d(3)*rcuti3 )
664 >             duduz_j(1) = duduz_j(1) - sw*pref*( ri2*xhat - d(1)*rcuti3 )
665 >             duduz_j(2) = duduz_j(2) - sw*pref*( ri2*yhat - d(2)*rcuti3 )
666 >             duduz_j(3) = duduz_j(3) - sw*pref*( ri2*zhat - d(3)*rcuti3 )
667  
668            else
669               if (j_is_SplitDipole) then
# Line 660 | Line 678 | contains
678               ri2 = ri * ri
679               ri3 = ri2 * ri
680               sc2 = scale * scale
681 <            
681 >
682 >             pref = pre12 * q_i * mu_j
683               vterm = - pref * ct_j * ri2 * scale
684 <             vpair = vpair + swi * vterm
685 <             epot = epot + vterm
684 >             vpair = vpair + vterm
685 >             epot = epot + sw*vterm
686              
687               !! this has a + sign in the () because the rij vector is
688               !! r_j - r_i and the charge-dipole potential takes the origin
689               !! as the point dipole, which is atom j in this case.
690              
691 <             dudx = dudx - pref * ri3 * ( uz_j(1) - 3.0d0*ct_j*xhat*sc2)
692 <             dudy = dudy - pref * ri3 * ( uz_j(2) - 3.0d0*ct_j*yhat*sc2)
693 <             dudz = dudz - pref * ri3 * ( uz_j(3) - 3.0d0*ct_j*zhat*sc2)
691 >             dudx = dudx - sw*pref * ri3 * ( uz_j(1) - 3.0d0*ct_j*xhat*sc2)
692 >             dudy = dudy - sw*pref * ri3 * ( uz_j(2) - 3.0d0*ct_j*yhat*sc2)
693 >             dudz = dudz - sw*pref * ri3 * ( uz_j(3) - 3.0d0*ct_j*zhat*sc2)
694              
695 <             duduz_j(1) = duduz_j(1) - pref * ri2 * xhat * scale
696 <             duduz_j(2) = duduz_j(2) - pref * ri2 * yhat * scale
697 <             duduz_j(3) = duduz_j(3) - pref * ri2 * zhat * scale
695 >             duduz_j(1) = duduz_j(1) - sw*pref * ri2 * xhat * scale
696 >             duduz_j(2) = duduz_j(2) - sw*pref * ri2 * yhat * scale
697 >             duduz_j(3) = duduz_j(3) - sw*pref * ri2 * zhat * scale
698  
699            endif
700         endif
# Line 688 | Line 707 | contains
707            cy2 = cy_j * cy_j
708            cz2 = cz_j * cz_j
709  
710 <
711 <          pref =  sw * pre14 * q_i / 3.0_dp
693 <
694 <          if (corrMethod .eq. 1) then
710 >          if (summationMethod .eq. UNDAMPED_WOLF) then
711 >             pref =  pre14 * q_i / 3.0_dp
712               vterm1 = pref * ri3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
713                    qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
714                    qzz_j * (3.0_dp*cz2 - 1.0_dp) )
715               vterm2 = pref * rcuti3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
716                    qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
717                    qzz_j * (3.0_dp*cz2 - 1.0_dp) )
718 <             vpair = vpair + swi*( vterm1 - vterm2 )
719 <             epot = epot + ( vterm1 - vterm2 )
718 >             vpair = vpair + ( vterm1 - vterm2 )
719 >             epot = epot + sw*( vterm1 - vterm2 )
720              
721               dudx = dudx - (5.0_dp * &
722 <                  (vterm1*riji*xhat - vterm2*rcuti2*d(1))) + pref * ( &
722 >                  (vterm1*riji*xhat - vterm2*rcuti2*d(1))) + sw*pref * ( &
723                    (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(1)) - &
724                    qxx_j*2.0_dp*(xhat - rcuti*d(1))) + &
725                    (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(1)) - &
# Line 710 | Line 727 | contains
727                    (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(1)) - &
728                    qzz_j*2.0_dp*(xhat - rcuti*d(1))) )
729               dudy = dudy - (5.0_dp * &
730 <                  (vterm1*riji*yhat - vterm2*rcuti2*d(2))) + pref * ( &
730 >                  (vterm1*riji*yhat - vterm2*rcuti2*d(2))) + sw*pref * ( &
731                    (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(2)) - &
732                    qxx_j*2.0_dp*(yhat - rcuti*d(2))) + &
733                    (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(2)) - &
# Line 718 | Line 735 | contains
735                    (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(2)) - &
736                    qzz_j*2.0_dp*(yhat - rcuti*d(2))) )
737               dudz = dudz - (5.0_dp * &
738 <                  (vterm1*riji*zhat - vterm2*rcuti2*d(3))) + pref * ( &
738 >                  (vterm1*riji*zhat - vterm2*rcuti2*d(3))) + sw*pref * ( &
739                    (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(3)) - &
740                    qxx_j*2.0_dp*(zhat - rcuti*d(3))) + &
741                    (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(3)) - &
# Line 726 | Line 743 | contains
743                    (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(3)) - &
744                    qzz_j*2.0_dp*(zhat - rcuti*d(3))) )
745              
746 <             dudux_j(1) = dudux_j(1) + pref * (ri3*(qxx_j*6.0_dp*cx_j*xhat) - &
746 >             dudux_j(1) = dudux_j(1) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*xhat) -&
747                    rcuti4*(qxx_j*6.0_dp*cx_j*d(1)))
748 <             dudux_j(2) = dudux_j(2) + pref * (ri3*(qxx_j*6.0_dp*cx_j*yhat) - &
748 >             dudux_j(2) = dudux_j(2) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*yhat) -&
749                    rcuti4*(qxx_j*6.0_dp*cx_j*d(2)))
750 <             dudux_j(3) = dudux_j(3) + pref * (ri3*(qxx_j*6.0_dp*cx_j*zhat) - &
750 >             dudux_j(3) = dudux_j(3) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*zhat) -&
751                    rcuti4*(qxx_j*6.0_dp*cx_j*d(3)))
752              
753 <             duduy_j(1) = duduy_j(1) + pref * (ri3*(qyy_j*6.0_dp*cy_j*xhat) - &
753 >             duduy_j(1) = duduy_j(1) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*xhat) -&
754                    rcuti4*(qyy_j*6.0_dp*cx_j*d(1)))
755 <             duduy_j(2) = duduy_j(2) + pref * (ri3*(qyy_j*6.0_dp*cy_j*yhat) - &
755 >             duduy_j(2) = duduy_j(2) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*yhat) -&
756                    rcuti4*(qyy_j*6.0_dp*cx_j*d(2)))
757 <             duduy_j(3) = duduy_j(3) + pref * (ri3*(qyy_j*6.0_dp*cy_j*zhat) - &
757 >             duduy_j(3) = duduy_j(3) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*zhat) -&
758                    rcuti4*(qyy_j*6.0_dp*cx_j*d(3)))
759              
760 <             duduz_j(1) = duduz_j(1) + pref * (ri3*(qzz_j*6.0_dp*cz_j*xhat) - &
760 >             duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*xhat) -&
761                    rcuti4*(qzz_j*6.0_dp*cx_j*d(1)))
762 <             duduz_j(2) = duduz_j(2) + pref * (ri3*(qzz_j*6.0_dp*cz_j*yhat) - &
762 >             duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*yhat) -&
763                    rcuti4*(qzz_j*6.0_dp*cx_j*d(2)))
764 <             duduz_j(3) = duduz_j(3) + pref * (ri3*(qzz_j*6.0_dp*cz_j*zhat) - &
764 >             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*zhat) -&
765                    rcuti4*(qzz_j*6.0_dp*cx_j*d(3)))
766          
767            else
768 +             pref =  pre14 * q_i / 3.0_dp
769               vterm = pref * ri3 * (qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
770                    qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
771                    qzz_j * (3.0_dp*cz2 - 1.0_dp))
772 <             vpair = vpair + swi * vterm
773 <             epot = epot + vterm
772 >             vpair = vpair + vterm
773 >             epot = epot + sw*vterm
774              
775 <             dudx = dudx - 5.0_dp*vterm*riji*xhat + pref * ri4 * ( &
775 >             dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref * ri4 * ( &
776                    qxx_j*(6.0_dp*cx_j*ux_j(1) - 2.0_dp*xhat) + &
777                    qyy_j*(6.0_dp*cy_j*uy_j(1) - 2.0_dp*xhat) + &
778                    qzz_j*(6.0_dp*cz_j*uz_j(1) - 2.0_dp*xhat) )
779 <             dudy = dudy - 5.0_dp*vterm*riji*yhat + pref * ri4 * ( &
779 >             dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref * ri4 * ( &
780                    qxx_j*(6.0_dp*cx_j*ux_j(2) - 2.0_dp*yhat) + &
781                    qyy_j*(6.0_dp*cy_j*uy_j(2) - 2.0_dp*yhat) + &
782                    qzz_j*(6.0_dp*cz_j*uz_j(2) - 2.0_dp*yhat) )
783 <             dudz = dudz - 5.0_dp*vterm*riji*zhat + pref * ri4 * ( &
783 >             dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref * ri4 * ( &
784                    qxx_j*(6.0_dp*cx_j*ux_j(3) - 2.0_dp*zhat) + &
785                    qyy_j*(6.0_dp*cy_j*uy_j(3) - 2.0_dp*zhat) + &
786                    qzz_j*(6.0_dp*cz_j*uz_j(3) - 2.0_dp*zhat) )
787              
788 <             dudux_j(1) = dudux_j(1) + pref * ri3*(qxx_j*6.0_dp*cx_j*xhat)
789 <             dudux_j(2) = dudux_j(2) + pref * ri3*(qxx_j*6.0_dp*cx_j*yhat)
790 <             dudux_j(3) = dudux_j(3) + pref * ri3*(qxx_j*6.0_dp*cx_j*zhat)
788 >             dudux_j(1) = dudux_j(1) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*xhat)
789 >             dudux_j(2) = dudux_j(2) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*yhat)
790 >             dudux_j(3) = dudux_j(3) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*zhat)
791              
792 <             duduy_j(1) = duduy_j(1) + pref * ri3*(qyy_j*6.0_dp*cy_j*xhat)
793 <             duduy_j(2) = duduy_j(2) + pref * ri3*(qyy_j*6.0_dp*cy_j*yhat)
794 <             duduy_j(3) = duduy_j(3) + pref * ri3*(qyy_j*6.0_dp*cy_j*zhat)
792 >             duduy_j(1) = duduy_j(1) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*xhat)
793 >             duduy_j(2) = duduy_j(2) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*yhat)
794 >             duduy_j(3) = duduy_j(3) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*zhat)
795              
796 <             duduz_j(1) = duduz_j(1) + pref * ri3*(qzz_j*6.0_dp*cz_j*xhat)
797 <             duduz_j(2) = duduz_j(2) + pref * ri3*(qzz_j*6.0_dp*cz_j*yhat)
798 <             duduz_j(3) = duduz_j(3) + pref * ri3*(qzz_j*6.0_dp*cz_j*zhat)
796 >             duduz_j(1) = duduz_j(1) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*xhat)
797 >             duduz_j(2) = duduz_j(2) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*yhat)
798 >             duduz_j(3) = duduz_j(3) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*zhat)
799            
800            endif
801         endif
# Line 786 | Line 804 | contains
804      if (i_is_Dipole) then
805  
806         if (j_is_Charge) then
807 <
808 <          pref = sw * pre12 * q_j * mu_i
809 <
810 <          if (corrMethod .eq. 1) then
807 >          
808 >          pref = pre12 * q_j * mu_i
809 >          
810 >          if (summationMethod .eq. UNDAMPED_WOLF) then
811               ri2 = riji * riji
812               ri3 = ri2 * riji
813  
814 +             pref = pre12 * q_j * mu_i
815               vterm = pref * ct_i * (ri2 - rcuti2)
816 <             vpair = vpair + swi * vterm
817 <             epot = epot + vterm
816 >             vpair = vpair + vterm
817 >             epot = epot + sw*vterm
818              
819               !! this has a + sign in the () because the rij vector is
820               !! r_j - r_i and the charge-dipole potential takes the origin
821               !! as the point dipole, which is atom j in this case.
822              
823 <             dudx = dudx + pref * ( ri3*( uz_i(1) - 3.0d0*ct_i*xhat) &
823 >             dudx = dudx + sw*pref * ( ri3*( uz_i(1) - 3.0d0*ct_i*xhat) &
824                    - rcuti3*( uz_i(1) - 3.0d0*ct_i*d(1)*rcuti ) )
825 <             dudy = dudy + pref * ( ri3*( uz_i(2) - 3.0d0*ct_i*yhat) &
825 >             dudy = dudy + sw*pref * ( ri3*( uz_i(2) - 3.0d0*ct_i*yhat) &
826                    - rcuti3*( uz_i(2) - 3.0d0*ct_i*d(2)*rcuti ) )
827 <             dudz = dudz + pref * ( ri3*( uz_i(3) - 3.0d0*ct_i*zhat) &
827 >             dudz = dudz + sw*pref * ( ri3*( uz_i(3) - 3.0d0*ct_i*zhat) &
828                    - rcuti3*( uz_i(3) - 3.0d0*ct_i*d(3)*rcuti ) )
829              
830 <             duduz_i(1) = duduz_i(1) - pref*( ri2*xhat - d(1)*rcuti3 )
831 <             duduz_i(2) = duduz_i(2) - pref*( ri2*yhat - d(2)*rcuti3 )
832 <             duduz_i(3) = duduz_i(3) - pref*( ri2*zhat - d(3)*rcuti3 )
830 >             duduz_i(1) = duduz_i(1) - sw*pref*( ri2*xhat - d(1)*rcuti3 )
831 >             duduz_i(2) = duduz_i(2) - sw*pref*( ri2*yhat - d(2)*rcuti3 )
832 >             duduz_i(3) = duduz_i(3) - sw*pref*( ri2*zhat - d(3)*rcuti3 )
833  
834            else
835               if (i_is_SplitDipole) then
# Line 825 | Line 844 | contains
844               ri2 = ri * ri
845               ri3 = ri2 * ri
846               sc2 = scale * scale
847 <            
847 >
848 >             pref = pre12 * q_j * mu_i
849               vterm = pref * ct_i * ri2 * scale
850 <             vpair = vpair + swi * vterm
851 <             epot = epot + vterm
850 >             vpair = vpair + vterm
851 >             epot = epot + sw*vterm
852              
853 <             dudx = dudx + pref * ri3 * ( uz_i(1) - 3.0d0 * ct_i * xhat*sc2)
854 <             dudy = dudy + pref * ri3 * ( uz_i(2) - 3.0d0 * ct_i * yhat*sc2)
855 <             dudz = dudz + pref * ri3 * ( uz_i(3) - 3.0d0 * ct_i * zhat*sc2)
853 >             dudx = dudx + sw*pref * ri3 * ( uz_i(1) - 3.0d0 * ct_i * xhat*sc2)
854 >             dudy = dudy + sw*pref * ri3 * ( uz_i(2) - 3.0d0 * ct_i * yhat*sc2)
855 >             dudz = dudz + sw*pref * ri3 * ( uz_i(3) - 3.0d0 * ct_i * zhat*sc2)
856              
857 <             duduz_i(1) = duduz_i(1) + pref * ri2 * xhat * scale
858 <             duduz_i(2) = duduz_i(2) + pref * ri2 * yhat * scale
859 <             duduz_i(3) = duduz_i(3) + pref * ri2 * zhat * scale
857 >             duduz_i(1) = duduz_i(1) + sw*pref * ri2 * xhat * scale
858 >             duduz_i(2) = duduz_i(2) + sw*pref * ri2 * yhat * scale
859 >             duduz_i(3) = duduz_i(3) + sw*pref * ri2 * zhat * scale
860            endif
861         endif
862 <
862 >      
863         if (j_is_Dipole) then
864  
865 <          pref = sw * pre22 * mu_i * mu_j
846 <
847 <          if (corrMethod .eq. 1) then
865 >          if (summationMethod .eq. UNDAMPED_WOLF) then
866               ri2 = riji * riji
867               ri3 = ri2 * riji
868               ri4 = ri2 * ri2
869  
870 +             pref = pre22 * mu_i * mu_j
871               vterm = pref * (ri3 - rcuti3) * (ct_ij - 3.0d0 * ct_i * ct_j)
872 <             vpair = vpair + swi * vterm
873 <             epot = epot + vterm
872 >             vpair = vpair + vterm
873 >             epot = epot + sw*vterm
874              
875               a1 = 5.0d0 * ct_i * ct_j - ct_ij
876              
877 <             dudx = dudx + pref*3.0d0*ri4 &
878 <                  *(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1)) - &
879 <                  pref*3.0d0*rcuti4*(a1*rcuti*d(1)-ct_i*uz_j(1)-ct_j*uz_i(1))
880 <             dudy = dudy + pref*3.0d0*ri4 &
881 <                  *(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2)) - &
882 <                  pref*3.0d0*rcuti4*(a1*rcuti*d(2)-ct_i*uz_j(2)-ct_j*uz_i(2))
883 <             dudz = dudz + pref*3.0d0*ri4 &
884 <                  *(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3)) - &
885 <                  pref*3.0d0*rcuti4*(a1*rcuti*d(3)-ct_i*uz_j(3)-ct_j*uz_i(3))
877 >             dudx = dudx + sw*pref*3.0d0*ri4 &
878 >                             * (a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1)) &
879 >                         - sw*pref*3.0d0*rcuti4 &
880 >                             * (a1*rcuti*d(1)-ct_i*uz_j(1)-ct_j*uz_i(1))
881 >             dudy = dudy + sw*pref*3.0d0*ri4 &
882 >                             * (a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2)) &
883 >                         - sw*pref*3.0d0*rcuti4 &
884 >                             * (a1*rcuti*d(2)-ct_i*uz_j(2)-ct_j*uz_i(2))
885 >             dudz = dudz + sw*pref*3.0d0*ri4 &
886 >                             * (a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3)) &
887 >                         - sw*pref*3.0d0*rcuti4 &
888 >                             * (a1*rcuti*d(3)-ct_i*uz_j(3)-ct_j*uz_i(3))
889              
890 <             duduz_i(1) = duduz_i(1) + pref*(ri3*(uz_j(1) - 3.0d0*ct_j*xhat) &
890 >             duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(uz_j(1)-3.0d0*ct_j*xhat) &
891                    - rcuti3*(uz_j(1) - 3.0d0*ct_j*d(1)*rcuti))
892 <             duduz_i(2) = duduz_i(2) + pref*(ri3*(uz_j(2) - 3.0d0*ct_j*yhat) &
892 >             duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(uz_j(2)-3.0d0*ct_j*yhat) &
893                    - rcuti3*(uz_j(2) - 3.0d0*ct_j*d(2)*rcuti))
894 <             duduz_i(3) = duduz_i(3) + pref*(ri3*(uz_j(3) - 3.0d0*ct_j*zhat) &
894 >             duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(uz_j(3)-3.0d0*ct_j*zhat) &
895                    - rcuti3*(uz_j(3) - 3.0d0*ct_j*d(3)*rcuti))
896 <             duduz_j(1) = duduz_j(1) + pref*(ri3*(uz_i(1) - 3.0d0*ct_i*xhat) &
896 >             duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(uz_i(1)-3.0d0*ct_i*xhat) &
897                    - rcuti3*(uz_i(1) - 3.0d0*ct_i*d(1)*rcuti))
898 <             duduz_j(2) = duduz_j(2) + pref*(ri3*(uz_i(2) - 3.0d0*ct_i*yhat) &
898 >             duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(uz_i(2)-3.0d0*ct_i*yhat) &
899                    - rcuti3*(uz_i(2) - 3.0d0*ct_i*d(2)*rcuti))
900 <             duduz_j(3) = duduz_j(3) + pref*(ri3*(uz_i(3) - 3.0d0*ct_i*zhat) &
900 >             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(uz_i(3)-3.0d0*ct_i*zhat) &
901                    - rcuti3*(uz_i(3) - 3.0d0*ct_i*d(3)*rcuti))
902 +
903            else
881 –            
904               if (i_is_SplitDipole) then
905                  if (j_is_SplitDipole) then
906                     BigR = sqrt(r2 + 0.25_dp * d_i * d_i + 0.25_dp * d_j * d_j)
# Line 905 | Line 927 | contains
927               ri4 = ri2 * ri2
928               sc2 = scale * scale
929              
930 +             pref = pre22 * mu_i * mu_j
931               vterm = pref * ri3 * (ct_ij - 3.0d0 * ct_i * ct_j * sc2)
932 <             vpair = vpair + swi * vterm
933 <             epot = epot + vterm
932 >             vpair = vpair + vterm
933 >             epot = epot + sw*vterm
934              
935               a1 = 5.0d0 * ct_i * ct_j * sc2 - ct_ij
936              
937 <             dudx = dudx + pref*3.0d0*ri4*scale &
938 <                  *(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
939 <             dudy = dudy + pref*3.0d0*ri4*scale &
940 <                  *(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
941 <             dudz = dudz + pref*3.0d0*ri4*scale &
942 <                  *(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
937 >             dudx = dudx + sw*pref*3.0d0*ri4*scale &
938 >                             *(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
939 >             dudy = dudy + sw*pref*3.0d0*ri4*scale &
940 >                             *(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
941 >             dudz = dudz + sw*pref*3.0d0*ri4*scale &
942 >                             *(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
943              
944 <             duduz_i(1) = duduz_i(1) + pref*ri3 &
945 <                  *(uz_j(1) - 3.0d0*ct_j*xhat*sc2)
946 <             duduz_i(2) = duduz_i(2) + pref*ri3 &
947 <                  *(uz_j(2) - 3.0d0*ct_j*yhat*sc2)
948 <             duduz_i(3) = duduz_i(3) + pref*ri3 &
949 <                  *(uz_j(3) - 3.0d0*ct_j*zhat*sc2)
944 >             duduz_i(1) = duduz_i(1) + sw*pref*ri3 &
945 >                                         *(uz_j(1) - 3.0d0*ct_j*xhat*sc2)
946 >             duduz_i(2) = duduz_i(2) + sw*pref*ri3 &
947 >                                         *(uz_j(2) - 3.0d0*ct_j*yhat*sc2)
948 >             duduz_i(3) = duduz_i(3) + sw*pref*ri3 &
949 >                                         *(uz_j(3) - 3.0d0*ct_j*zhat*sc2)
950              
951 <             duduz_j(1) = duduz_j(1) + pref*ri3 &
952 <                  *(uz_i(1) - 3.0d0*ct_i*xhat*sc2)
953 <             duduz_j(2) = duduz_j(2) + pref*ri3 &
954 <                  *(uz_i(2) - 3.0d0*ct_i*yhat*sc2)
955 <             duduz_j(3) = duduz_j(3) + pref*ri3 &
956 <                  *(uz_i(3) - 3.0d0*ct_i*zhat*sc2)
951 >             duduz_j(1) = duduz_j(1) + sw*pref*ri3 &
952 >                                         *(uz_i(1) - 3.0d0*ct_i*xhat*sc2)
953 >             duduz_j(2) = duduz_j(2) + sw*pref*ri3 &
954 >                                         *(uz_i(2) - 3.0d0*ct_i*yhat*sc2)
955 >             duduz_j(3) = duduz_j(3) + sw*pref*ri3 &
956 >                                         *(uz_i(3) - 3.0d0*ct_i*zhat*sc2)
957            endif
958         endif
959      endif
# Line 945 | Line 968 | contains
968            cy2 = cy_i * cy_i
969            cz2 = cz_i * cz_i
970  
971 <          pref = sw * pre14 * q_j / 3.0_dp
972 <
950 <          if (corrMethod .eq. 1) then
971 >          if (summationMethod .eq. UNDAMPED_WOLF) then
972 >             pref = pre14 * q_j / 3.0_dp
973               vterm1 = pref * ri3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
974                    qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
975                    qzz_i * (3.0_dp*cz2 - 1.0_dp) )
976               vterm2 = pref * rcuti3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
977                    qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
978                    qzz_i * (3.0_dp*cz2 - 1.0_dp) )
979 <             vpair = vpair + swi * ( vterm1 - vterm2 )
980 <             epot = epot + ( vterm1 - vterm2 )
979 >             vpair = vpair + ( vterm1 - vterm2 )
980 >             epot = epot + sw*( vterm1 - vterm2 )
981              
982 <             dudx = dudx - (5.0_dp*(vterm1*riji*xhat - vterm2*rcuti2*d(1))) + &
983 <                  pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(1)) - &
982 >             dudx = dudx - sw*(5.0_dp*(vterm1*riji*xhat-vterm2*rcuti2*d(1))) +&
983 >                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(1)) - &
984                    qxx_i*2.0_dp*(xhat - rcuti*d(1))) + &
985                    (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(1)) - &
986                    qyy_i*2.0_dp*(xhat - rcuti*d(1))) + &
987                    (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(1)) - &
988                    qzz_i*2.0_dp*(xhat - rcuti*d(1))) )
989 <             dudy = dudy - (5.0_dp*(vterm1*riji*yhat - vterm2*rcuti2*d(2))) + &
990 <                  pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(2)) - &
989 >             dudy = dudy - sw*(5.0_dp*(vterm1*riji*yhat-vterm2*rcuti2*d(2))) +&
990 >                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(2)) - &
991                    qxx_i*2.0_dp*(yhat - rcuti*d(2))) + &
992                    (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(2)) - &
993                    qyy_i*2.0_dp*(yhat - rcuti*d(2))) + &
994                    (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(2)) - &
995                    qzz_i*2.0_dp*(yhat - rcuti*d(2))) )
996 <             dudz = dudz - (5.0_dp*(vterm1*riji*zhat - vterm2*rcuti2*d(3))) + &
997 <                  pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(3)) - &
996 >             dudz = dudz - sw*(5.0_dp*(vterm1*riji*zhat-vterm2*rcuti2*d(3))) +&
997 >                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(3)) - &
998                    qxx_i*2.0_dp*(zhat - rcuti*d(3))) + &
999                    (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(3)) - &
1000                    qyy_i*2.0_dp*(zhat - rcuti*d(3))) + &
1001                    (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(3)) - &
1002                    qzz_i*2.0_dp*(zhat - rcuti*d(3))) )
1003              
1004 <             dudux_i(1) = dudux_i(1) + pref * (ri3*(qxx_i*6.0_dp*cx_i*xhat) - &
1004 >             dudux_i(1) = dudux_i(1) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*xhat) -&
1005                    rcuti4*(qxx_i*6.0_dp*cx_i*d(1)))
1006 <             dudux_i(2) = dudux_i(2) + pref * (ri3*(qxx_i*6.0_dp*cx_i*yhat) - &
1006 >             dudux_i(2) = dudux_i(2) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*yhat) -&
1007                    rcuti4*(qxx_i*6.0_dp*cx_i*d(2)))
1008 <             dudux_i(3) = dudux_i(3) + pref * (ri3*(qxx_i*6.0_dp*cx_i*zhat) - &
1008 >             dudux_i(3) = dudux_i(3) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*zhat) -&
1009                    rcuti4*(qxx_i*6.0_dp*cx_i*d(3)))
1010              
1011 <             duduy_i(1) = duduy_i(1) + pref * (ri3*(qyy_i*6.0_dp*cy_i*xhat) - &
1011 >             duduy_i(1) = duduy_i(1) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*xhat) -&
1012                    rcuti4*(qyy_i*6.0_dp*cx_i*d(1)))
1013 <             duduy_i(2) = duduy_i(2) + pref * (ri3*(qyy_i*6.0_dp*cy_i*yhat) - &
1013 >             duduy_i(2) = duduy_i(2) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*yhat) -&
1014                    rcuti4*(qyy_i*6.0_dp*cx_i*d(2)))
1015 <             duduy_i(3) = duduy_i(3) + pref * (ri3*(qyy_i*6.0_dp*cy_i*zhat) - &
1015 >             duduy_i(3) = duduy_i(3) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*zhat) -&
1016                    rcuti4*(qyy_i*6.0_dp*cx_i*d(3)))
1017              
1018 <             duduz_i(1) = duduz_i(1) + pref * (ri3*(qzz_i*6.0_dp*cz_i*xhat) - &
1018 >             duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*xhat) -&
1019                    rcuti4*(qzz_i*6.0_dp*cx_i*d(1)))
1020 <             duduz_i(2) = duduz_i(2) + pref * (ri3*(qzz_i*6.0_dp*cz_i*yhat) - &
1020 >             duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*yhat) -&
1021                    rcuti4*(qzz_i*6.0_dp*cx_i*d(2)))
1022 <             duduz_i(3) = duduz_i(3) + pref * (ri3*(qzz_i*6.0_dp*cz_i*zhat) - &
1022 >             duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*zhat) -&
1023                    rcuti4*(qzz_i*6.0_dp*cx_i*d(3)))
1024  
1025            else
1026 +             pref = pre14 * q_j / 3.0_dp
1027               vterm = pref * ri3 * (qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1028                    qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1029                    qzz_i * (3.0_dp*cz2 - 1.0_dp))
1030 <             vpair = vpair + swi * vterm
1031 <             epot = epot + vterm
1030 >             vpair = vpair + vterm
1031 >             epot = epot + sw*vterm
1032              
1033 <             dudx = dudx - 5.0_dp*vterm*riji*xhat + pref * ri4 * ( &
1033 >             dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref*ri4 * ( &
1034                    qxx_i*(6.0_dp*cx_i*ux_i(1) - 2.0_dp*xhat) + &
1035                    qyy_i*(6.0_dp*cy_i*uy_i(1) - 2.0_dp*xhat) + &
1036                    qzz_i*(6.0_dp*cz_i*uz_i(1) - 2.0_dp*xhat) )
1037 <             dudy = dudy - 5.0_dp*vterm*riji*yhat + pref * ri4 * ( &
1037 >             dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref*ri4 * ( &
1038                    qxx_i*(6.0_dp*cx_i*ux_i(2) - 2.0_dp*yhat) + &
1039                    qyy_i*(6.0_dp*cy_i*uy_i(2) - 2.0_dp*yhat) + &
1040                    qzz_i*(6.0_dp*cz_i*uz_i(2) - 2.0_dp*yhat) )
1041 <             dudz = dudz - 5.0_dp*vterm*riji*zhat + pref * ri4 * ( &
1041 >             dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref*ri4 * ( &
1042                    qxx_i*(6.0_dp*cx_i*ux_i(3) - 2.0_dp*zhat) + &
1043                    qyy_i*(6.0_dp*cy_i*uy_i(3) - 2.0_dp*zhat) + &
1044                    qzz_i*(6.0_dp*cz_i*uz_i(3) - 2.0_dp*zhat) )
1045              
1046 <             dudux_i(1) = dudux_i(1) + pref * ri3*(qxx_i*6.0_dp*cx_i*xhat)
1047 <             dudux_i(2) = dudux_i(2) + pref * ri3*(qxx_i*6.0_dp*cx_i*yhat)
1048 <             dudux_i(3) = dudux_i(3) + pref * ri3*(qxx_i*6.0_dp*cx_i*zhat)
1046 >             dudux_i(1) = dudux_i(1) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*xhat)
1047 >             dudux_i(2) = dudux_i(2) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*yhat)
1048 >             dudux_i(3) = dudux_i(3) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*zhat)
1049              
1050 <             duduy_i(1) = duduy_i(1) + pref * ri3*(qyy_i*6.0_dp*cy_i*xhat)
1051 <             duduy_i(2) = duduy_i(2) + pref * ri3*(qyy_i*6.0_dp*cy_i*yhat)
1052 <             duduy_i(3) = duduy_i(3) + pref * ri3*(qyy_i*6.0_dp*cy_i*zhat)
1050 >             duduy_i(1) = duduy_i(1) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*xhat)
1051 >             duduy_i(2) = duduy_i(2) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*yhat)
1052 >             duduy_i(3) = duduy_i(3) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*zhat)
1053              
1054 <             duduz_i(1) = duduz_i(1) + pref * ri3*(qzz_i*6.0_dp*cz_i*xhat)
1055 <             duduz_i(2) = duduz_i(2) + pref * ri3*(qzz_i*6.0_dp*cz_i*yhat)
1056 <             duduz_i(3) = duduz_i(3) + pref * ri3*(qzz_i*6.0_dp*cz_i*zhat)
1054 >             duduz_i(1) = duduz_i(1) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*xhat)
1055 >             duduz_i(2) = duduz_i(2) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*yhat)
1056 >             duduz_i(3) = duduz_i(3) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*zhat)
1057            endif
1058         endif
1059      endif

Diff Legend

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