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 2302 by chrisfen, Fri Sep 16 16:07:39 2005 UTC vs.
Revision 2343 by chrisfen, Tue Oct 4 19:33:22 2005 UTC

# Line 83 | 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 <
87 <
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 >  real(kind=dp), save :: alphaPi = 0.0_dp
91 >  real(kind=dp), save :: invRootPi = 0.0_dp
92 >  
93 > #ifdef __IFC
94 > ! error function for ifc version > 7.
95 >  double precision, external :: derfc
96 > #endif
97 >  
98    public :: setElectrostaticSummationMethod
99    public :: setElectrostaticCutoffRadius
100    public :: setDampedWolfAlpha
# Line 118 | Line 128 | contains
128   contains
129  
130    subroutine setElectrostaticSummationMethod(the_ESM)
121 –
131      integer, intent(in) :: the_ESM    
132  
133      if ((the_ESM .le. 0) .or. (the_ESM .gt. REACTION_FIELD)) then
134         call handleError("setElectrostaticSummationMethod", "Unsupported Summation Method")
135      endif
136  
137 +    summationMethod = the_ESM
138 +
139    end subroutine setElectrostaticSummationMethod
140  
141    subroutine setElectrostaticCutoffRadius(thisRcut)
# Line 357 | Line 368 | contains
368  
369    subroutine checkSummationMethod()
370  
371 +    if (.not.haveDefaultCutoff) then
372 +       call handleError("checkSummationMethod", "no Default Cutoff set!")
373 +    endif
374 +
375 +    rcuti = 1.0d0 / defaultCutoff
376 +    rcuti2 = rcuti*rcuti
377 +    rcuti3 = rcuti2*rcuti
378 +    rcuti4 = rcuti2*rcuti2
379 +
380      if (summationMethod .eq. DAMPED_WOLF) then
381         if (.not.haveDWAconstants) then
382            
# Line 369 | Line 389 | contains
389            endif
390  
391            constEXP = exp(-dampingAlpha*dampingAlpha*defaultCutoff*defaultCutoff)
392 <          constERFC = erfc(dampingAlpha*defaultCutoff)
393 <          
392 >          constERFC = derfc(dampingAlpha*defaultCutoff)
393 >          invRootPi = 0.56418958354775628695d0
394 >          alphaPi = 2*dampingAlpha*invRootPi
395 >  
396            haveDWAconstants = .true.
397         endif
398      endif
# Line 387 | Line 409 | contains
409  
410  
411    subroutine doElectrostaticPair(atom1, atom2, d, rij, r2, sw, &
412 <       vpair, fpair, pot, eFrame, f, t, do_pot, corrMethod, rcuti)
412 >       vpair, fpair, pot, eFrame, f, t, do_pot)
413  
414      logical, intent(in) :: do_pot
415  
416      integer, intent(in) :: atom1, atom2
417      integer :: localError
396 –    integer, intent(in) :: corrMethod
418  
419 <    real(kind=dp), intent(in) :: rij, r2, sw, rcuti
419 >    real(kind=dp), intent(in) :: rij, r2, sw
420      real(kind=dp), intent(in), dimension(3) :: d
421      real(kind=dp), intent(inout) :: vpair
422      real(kind=dp), intent(inout), dimension(3) :: fpair
423  
424 <    real( kind = dp ) :: pot, swi
424 >    real( kind = dp ) :: pot
425      real( kind = dp ), dimension(9,nLocal) :: eFrame
426      real( kind = dp ), dimension(3,nLocal) :: f
427      real( kind = dp ), dimension(3,nLocal) :: t
# Line 425 | Line 446 | contains
446      real (kind=dp) :: pref, vterm, epot, dudr, vterm1, vterm2
447      real (kind=dp) :: xhat, yhat, zhat
448      real (kind=dp) :: dudx, dudy, dudz
449 <    real (kind=dp) :: scale, sc2, bigR, switcher, dswitcher
450 <    real (kind=dp) :: rcuti2, rcuti3, rcuti4
449 >    real (kind=dp) :: scale, sc2, bigR
450 >    real (kind=dp) :: varERFC, varEXP
451 >    real (kind=dp) :: limScale
452  
453      if (.not.allocated(ElectrostaticMap)) then
454         call handleError("electrostatic", "no ElectrostaticMap was present before first call of do_electrostatic_pair!")
# Line 435 | Line 457 | contains
457  
458      if (.not.summationMethodChecked) then
459         call checkSummationMethod()
460 +      
461      endif
462  
463  
# Line 449 | Line 472 | contains
472      !! some variables we'll need independent of electrostatic type:
473  
474      riji = 1.0d0 / rij
475 <
475 >  
476      xhat = d(1) * riji
477      yhat = d(2) * riji
478      zhat = d(3) * riji
479  
457 –    rcuti2 = rcuti*rcuti
458 –    rcuti3 = rcuti2*rcuti
459 –    rcuti4 = rcuti2*rcuti2
460 –
461 –    swi = 1.0d0 / sw
462 –
480      !! logicals
481      i_is_Charge = ElectrostaticMap(me1)%is_Charge
482      i_is_Dipole = ElectrostaticMap(me1)%is_Dipole
# Line 578 | Line 595 | contains
595         cz_j = uz_j(1)*xhat + uz_j(2)*yhat + uz_j(3)*zhat
596      endif
597    
581 – !!$    switcher = 1.0d0
582 – !!$    dswitcher = 0.0d0
583 – !!$    ebalance = 0.0d0
584 – !!$    ! weaken the dipole interaction at close range for TAP water
585 – !!$    if (j_is_Tap .and. i_is_Tap) then
586 – !!$      call calc_switch(rij, mu_i, switcher, dswitcher)
587 – !!$    endif
588 –
598      epot = 0.0_dp
599      dudx = 0.0_dp
600      dudy = 0.0_dp
# Line 603 | Line 612 | contains
612  
613         if (j_is_Charge) then
614  
615 <          if (corrMethod .eq. 1) then
615 >          if (summationMethod .eq. UNDAMPED_WOLF) then
616 >
617               vterm = pre11 * q_i * q_j * (riji - rcuti)
618 +             vpair = vpair + vterm
619 +             epot = epot + sw*vterm
620 +            
621 +             dudr  = -sw*pre11*q_i*q_j * (riji*riji-rcuti2)*riji
622 +            
623 +             dudx = dudx + dudr * d(1)
624 +             dudy = dudy + dudr * d(2)
625 +             dudz = dudz + dudr * d(3)
626  
627 +          elseif (summationMethod .eq. DAMPED_WOLF) then
628 +
629 +             varERFC = derfc(dampingAlpha*rij)
630 +             varEXP = exp(-dampingAlpha*dampingAlpha*rij*rij)
631 +             vterm = pre11 * q_i * q_j * (varERFC*riji - constERFC*rcuti)
632               vpair = vpair + vterm
633 <             epot = epot + sw * vterm
633 >             epot = epot + sw*vterm
634              
635 <             dudr  = - sw * pre11 * q_i * q_j * (riji*riji*riji - rcuti2*rcuti)
635 >             dudr  = -sw*pre11*q_i*q_j * ( riji*((varERFC*riji*riji &
636 >                                                  + alphaPi*varEXP) &
637 >                                                 - (constERFC*rcuti2 &
638 >                                                    + alphaPi*constEXP)) )
639              
640               dudx = dudx + dudr * d(1)
641               dudy = dudy + dudr * d(2)
642               dudz = dudz + dudr * d(3)
643  
644            else
619 –             vterm = pre11 * q_i * q_j * riji
645  
646 +             vterm = pre11 * q_i * q_j * riji
647               vpair = vpair + vterm
648 <             epot = epot + sw * vterm
648 >             epot = epot + sw*vterm
649              
650               dudr  = - sw * vterm * riji
651              
# Line 633 | Line 659 | contains
659  
660         if (j_is_Dipole) then
661  
662 <          pref = sw * pre12 * q_i * mu_j
662 >          pref = pre12 * q_i * mu_j
663  
664 <          if (corrMethod .eq. 1) then
664 >          if (summationMethod .eq. UNDAMPED_WOLF) then
665               ri2 = riji * riji
666               ri3 = ri2 * riji
667  
668 +             pref = pre12 * q_i * mu_j
669               vterm = - pref * ct_j * (ri2 - rcuti2)
670 <             vpair = vpair + swi*vterm
671 <             epot = epot + vterm
670 >             vpair = vpair + vterm
671 >             epot = epot + sw*vterm
672              
673               !! this has a + sign in the () because the rij vector is
674               !! r_j - r_i and the charge-dipole potential takes the origin
675               !! as the point dipole, which is atom j in this case.
676              
677 <             dudx = dudx - pref * ( ri3*( uz_j(1) - 3.0d0*ct_j*xhat) &
677 >             dudx = dudx - sw*pref * ( ri3*( uz_j(1) - 3.0d0*ct_j*xhat) &
678                    - rcuti3*( uz_j(1) - 3.0d0*ct_j*d(1)*rcuti ) )
679 <             dudy = dudy - pref * ( ri3*( uz_j(2) - 3.0d0*ct_j*yhat) &
679 >             dudy = dudy - sw*pref * ( ri3*( uz_j(2) - 3.0d0*ct_j*yhat) &
680                    - rcuti3*( uz_j(2) - 3.0d0*ct_j*d(2)*rcuti ) )
681 <             dudz = dudz - pref * ( ri3*( uz_j(3) - 3.0d0*ct_j*zhat) &
681 >             dudz = dudz - sw*pref * ( ri3*( uz_j(3) - 3.0d0*ct_j*zhat) &
682                    - rcuti3*( uz_j(3) - 3.0d0*ct_j*d(3)*rcuti ) )
683              
684 <             duduz_j(1) = duduz_j(1) - pref*( ri2*xhat - d(1)*rcuti3 )
685 <             duduz_j(2) = duduz_j(2) - pref*( ri2*yhat - d(2)*rcuti3 )
686 <             duduz_j(3) = duduz_j(3) - pref*( ri2*zhat - d(3)*rcuti3 )
684 >             duduz_j(1) = duduz_j(1) - sw*pref*( ri2*xhat - d(1)*rcuti3 )
685 >             duduz_j(2) = duduz_j(2) - sw*pref*( ri2*yhat - d(2)*rcuti3 )
686 >             duduz_j(3) = duduz_j(3) - sw*pref*( ri2*zhat - d(3)*rcuti3 )
687  
688            else
689               if (j_is_SplitDipole) then
# Line 671 | Line 698 | contains
698               ri2 = ri * ri
699               ri3 = ri2 * ri
700               sc2 = scale * scale
701 <            
701 >
702 >             pref = pre12 * q_i * mu_j
703               vterm = - pref * ct_j * ri2 * scale
704 <             vpair = vpair + swi * vterm
705 <             epot = epot + vterm
704 >             vpair = vpair + vterm
705 >             epot = epot + sw*vterm
706              
707               !! this has a + sign in the () because the rij vector is
708               !! r_j - r_i and the charge-dipole potential takes the origin
709               !! as the point dipole, which is atom j in this case.
710              
711 <             dudx = dudx - pref * ri3 * ( uz_j(1) - 3.0d0*ct_j*xhat*sc2)
712 <             dudy = dudy - pref * ri3 * ( uz_j(2) - 3.0d0*ct_j*yhat*sc2)
713 <             dudz = dudz - pref * ri3 * ( uz_j(3) - 3.0d0*ct_j*zhat*sc2)
711 >             dudx = dudx - sw*pref * ri3 * ( uz_j(1) - 3.0d0*ct_j*xhat*sc2)
712 >             dudy = dudy - sw*pref * ri3 * ( uz_j(2) - 3.0d0*ct_j*yhat*sc2)
713 >             dudz = dudz - sw*pref * ri3 * ( uz_j(3) - 3.0d0*ct_j*zhat*sc2)
714              
715 <             duduz_j(1) = duduz_j(1) - pref * ri2 * xhat * scale
716 <             duduz_j(2) = duduz_j(2) - pref * ri2 * yhat * scale
717 <             duduz_j(3) = duduz_j(3) - pref * ri2 * zhat * scale
715 >             duduz_j(1) = duduz_j(1) - sw*pref * ri2 * xhat * scale
716 >             duduz_j(2) = duduz_j(2) - sw*pref * ri2 * yhat * scale
717 >             duduz_j(3) = duduz_j(3) - sw*pref * ri2 * zhat * scale
718  
719            endif
720         endif
# Line 699 | Line 727 | contains
727            cy2 = cy_j * cy_j
728            cz2 = cz_j * cz_j
729  
730 <
731 <          pref =  sw * pre14 * q_i / 3.0_dp
704 <
705 <          if (corrMethod .eq. 1) then
730 >          if (summationMethod .eq. UNDAMPED_WOLF) then
731 >             pref =  pre14 * q_i / 3.0_dp
732               vterm1 = pref * ri3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
733                    qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
734                    qzz_j * (3.0_dp*cz2 - 1.0_dp) )
735               vterm2 = pref * rcuti3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
736                    qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
737                    qzz_j * (3.0_dp*cz2 - 1.0_dp) )
738 <             vpair = vpair + swi*( vterm1 - vterm2 )
739 <             epot = epot + ( vterm1 - vterm2 )
738 >             vpair = vpair + ( vterm1 - vterm2 )
739 >             epot = epot + sw*( vterm1 - vterm2 )
740              
741               dudx = dudx - (5.0_dp * &
742 <                  (vterm1*riji*xhat - vterm2*rcuti2*d(1))) + pref * ( &
742 >                  (vterm1*riji*xhat - vterm2*rcuti2*d(1))) + sw*pref * ( &
743                    (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(1)) - &
744                    qxx_j*2.0_dp*(xhat - rcuti*d(1))) + &
745                    (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(1)) - &
# Line 721 | Line 747 | contains
747                    (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(1)) - &
748                    qzz_j*2.0_dp*(xhat - rcuti*d(1))) )
749               dudy = dudy - (5.0_dp * &
750 <                  (vterm1*riji*yhat - vterm2*rcuti2*d(2))) + pref * ( &
750 >                  (vterm1*riji*yhat - vterm2*rcuti2*d(2))) + sw*pref * ( &
751                    (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(2)) - &
752                    qxx_j*2.0_dp*(yhat - rcuti*d(2))) + &
753                    (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(2)) - &
# Line 729 | Line 755 | contains
755                    (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(2)) - &
756                    qzz_j*2.0_dp*(yhat - rcuti*d(2))) )
757               dudz = dudz - (5.0_dp * &
758 <                  (vterm1*riji*zhat - vterm2*rcuti2*d(3))) + pref * ( &
758 >                  (vterm1*riji*zhat - vterm2*rcuti2*d(3))) + sw*pref * ( &
759                    (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(3)) - &
760                    qxx_j*2.0_dp*(zhat - rcuti*d(3))) + &
761                    (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(3)) - &
# Line 737 | Line 763 | contains
763                    (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(3)) - &
764                    qzz_j*2.0_dp*(zhat - rcuti*d(3))) )
765              
766 <             dudux_j(1) = dudux_j(1) + pref * (ri3*(qxx_j*6.0_dp*cx_j*xhat) - &
766 >             dudux_j(1) = dudux_j(1) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*xhat) -&
767                    rcuti4*(qxx_j*6.0_dp*cx_j*d(1)))
768 <             dudux_j(2) = dudux_j(2) + pref * (ri3*(qxx_j*6.0_dp*cx_j*yhat) - &
768 >             dudux_j(2) = dudux_j(2) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*yhat) -&
769                    rcuti4*(qxx_j*6.0_dp*cx_j*d(2)))
770 <             dudux_j(3) = dudux_j(3) + pref * (ri3*(qxx_j*6.0_dp*cx_j*zhat) - &
770 >             dudux_j(3) = dudux_j(3) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*zhat) -&
771                    rcuti4*(qxx_j*6.0_dp*cx_j*d(3)))
772              
773 <             duduy_j(1) = duduy_j(1) + pref * (ri3*(qyy_j*6.0_dp*cy_j*xhat) - &
773 >             duduy_j(1) = duduy_j(1) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*xhat) -&
774                    rcuti4*(qyy_j*6.0_dp*cx_j*d(1)))
775 <             duduy_j(2) = duduy_j(2) + pref * (ri3*(qyy_j*6.0_dp*cy_j*yhat) - &
775 >             duduy_j(2) = duduy_j(2) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*yhat) -&
776                    rcuti4*(qyy_j*6.0_dp*cx_j*d(2)))
777 <             duduy_j(3) = duduy_j(3) + pref * (ri3*(qyy_j*6.0_dp*cy_j*zhat) - &
777 >             duduy_j(3) = duduy_j(3) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*zhat) -&
778                    rcuti4*(qyy_j*6.0_dp*cx_j*d(3)))
779              
780 <             duduz_j(1) = duduz_j(1) + pref * (ri3*(qzz_j*6.0_dp*cz_j*xhat) - &
780 >             duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*xhat) -&
781                    rcuti4*(qzz_j*6.0_dp*cx_j*d(1)))
782 <             duduz_j(2) = duduz_j(2) + pref * (ri3*(qzz_j*6.0_dp*cz_j*yhat) - &
782 >             duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*yhat) -&
783                    rcuti4*(qzz_j*6.0_dp*cx_j*d(2)))
784 <             duduz_j(3) = duduz_j(3) + pref * (ri3*(qzz_j*6.0_dp*cz_j*zhat) - &
784 >             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*zhat) -&
785                    rcuti4*(qzz_j*6.0_dp*cx_j*d(3)))
786          
787            else
788 +             pref =  pre14 * q_i / 3.0_dp
789               vterm = pref * ri3 * (qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
790                    qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
791                    qzz_j * (3.0_dp*cz2 - 1.0_dp))
792 <             vpair = vpair + swi * vterm
793 <             epot = epot + vterm
792 >             vpair = vpair + vterm
793 >             epot = epot + sw*vterm
794              
795 <             dudx = dudx - 5.0_dp*vterm*riji*xhat + pref * ri4 * ( &
795 >             dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref * ri4 * ( &
796                    qxx_j*(6.0_dp*cx_j*ux_j(1) - 2.0_dp*xhat) + &
797                    qyy_j*(6.0_dp*cy_j*uy_j(1) - 2.0_dp*xhat) + &
798                    qzz_j*(6.0_dp*cz_j*uz_j(1) - 2.0_dp*xhat) )
799 <             dudy = dudy - 5.0_dp*vterm*riji*yhat + pref * ri4 * ( &
799 >             dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref * ri4 * ( &
800                    qxx_j*(6.0_dp*cx_j*ux_j(2) - 2.0_dp*yhat) + &
801                    qyy_j*(6.0_dp*cy_j*uy_j(2) - 2.0_dp*yhat) + &
802                    qzz_j*(6.0_dp*cz_j*uz_j(2) - 2.0_dp*yhat) )
803 <             dudz = dudz - 5.0_dp*vterm*riji*zhat + pref * ri4 * ( &
803 >             dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref * ri4 * ( &
804                    qxx_j*(6.0_dp*cx_j*ux_j(3) - 2.0_dp*zhat) + &
805                    qyy_j*(6.0_dp*cy_j*uy_j(3) - 2.0_dp*zhat) + &
806                    qzz_j*(6.0_dp*cz_j*uz_j(3) - 2.0_dp*zhat) )
807              
808 <             dudux_j(1) = dudux_j(1) + pref * ri3*(qxx_j*6.0_dp*cx_j*xhat)
809 <             dudux_j(2) = dudux_j(2) + pref * ri3*(qxx_j*6.0_dp*cx_j*yhat)
810 <             dudux_j(3) = dudux_j(3) + pref * ri3*(qxx_j*6.0_dp*cx_j*zhat)
808 >             dudux_j(1) = dudux_j(1) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*xhat)
809 >             dudux_j(2) = dudux_j(2) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*yhat)
810 >             dudux_j(3) = dudux_j(3) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*zhat)
811              
812 <             duduy_j(1) = duduy_j(1) + pref * ri3*(qyy_j*6.0_dp*cy_j*xhat)
813 <             duduy_j(2) = duduy_j(2) + pref * ri3*(qyy_j*6.0_dp*cy_j*yhat)
814 <             duduy_j(3) = duduy_j(3) + pref * ri3*(qyy_j*6.0_dp*cy_j*zhat)
812 >             duduy_j(1) = duduy_j(1) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*xhat)
813 >             duduy_j(2) = duduy_j(2) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*yhat)
814 >             duduy_j(3) = duduy_j(3) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*zhat)
815              
816 <             duduz_j(1) = duduz_j(1) + pref * ri3*(qzz_j*6.0_dp*cz_j*xhat)
817 <             duduz_j(2) = duduz_j(2) + pref * ri3*(qzz_j*6.0_dp*cz_j*yhat)
818 <             duduz_j(3) = duduz_j(3) + pref * ri3*(qzz_j*6.0_dp*cz_j*zhat)
816 >             duduz_j(1) = duduz_j(1) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*xhat)
817 >             duduz_j(2) = duduz_j(2) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*yhat)
818 >             duduz_j(3) = duduz_j(3) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*zhat)
819            
820            endif
821         endif
# Line 797 | Line 824 | contains
824      if (i_is_Dipole) then
825  
826         if (j_is_Charge) then
827 <
828 <          pref = sw * pre12 * q_j * mu_i
829 <
830 <          if (corrMethod .eq. 1) then
827 >          
828 >          pref = pre12 * q_j * mu_i
829 >          
830 >          if (summationMethod .eq. UNDAMPED_WOLF) then
831               ri2 = riji * riji
832               ri3 = ri2 * riji
833  
834 +             pref = pre12 * q_j * mu_i
835               vterm = pref * ct_i * (ri2 - rcuti2)
836 <             vpair = vpair + swi * vterm
837 <             epot = epot + vterm
836 >             vpair = vpair + vterm
837 >             epot = epot + sw*vterm
838              
839               !! this has a + sign in the () because the rij vector is
840               !! r_j - r_i and the charge-dipole potential takes the origin
841               !! as the point dipole, which is atom j in this case.
842              
843 <             dudx = dudx + pref * ( ri3*( uz_i(1) - 3.0d0*ct_i*xhat) &
843 >             dudx = dudx + sw*pref * ( ri3*( uz_i(1) - 3.0d0*ct_i*xhat) &
844                    - rcuti3*( uz_i(1) - 3.0d0*ct_i*d(1)*rcuti ) )
845 <             dudy = dudy + pref * ( ri3*( uz_i(2) - 3.0d0*ct_i*yhat) &
845 >             dudy = dudy + sw*pref * ( ri3*( uz_i(2) - 3.0d0*ct_i*yhat) &
846                    - rcuti3*( uz_i(2) - 3.0d0*ct_i*d(2)*rcuti ) )
847 <             dudz = dudz + pref * ( ri3*( uz_i(3) - 3.0d0*ct_i*zhat) &
847 >             dudz = dudz + sw*pref * ( ri3*( uz_i(3) - 3.0d0*ct_i*zhat) &
848                    - rcuti3*( uz_i(3) - 3.0d0*ct_i*d(3)*rcuti ) )
849              
850 <             duduz_i(1) = duduz_i(1) - pref*( ri2*xhat - d(1)*rcuti3 )
851 <             duduz_i(2) = duduz_i(2) - pref*( ri2*yhat - d(2)*rcuti3 )
852 <             duduz_i(3) = duduz_i(3) - pref*( ri2*zhat - d(3)*rcuti3 )
850 >             duduz_i(1) = duduz_i(1) - sw*pref*( ri2*xhat - d(1)*rcuti3 )
851 >             duduz_i(2) = duduz_i(2) - sw*pref*( ri2*yhat - d(2)*rcuti3 )
852 >             duduz_i(3) = duduz_i(3) - sw*pref*( ri2*zhat - d(3)*rcuti3 )
853  
854            else
855               if (i_is_SplitDipole) then
# Line 836 | Line 864 | contains
864               ri2 = ri * ri
865               ri3 = ri2 * ri
866               sc2 = scale * scale
867 <            
867 >
868 >             pref = pre12 * q_j * mu_i
869               vterm = pref * ct_i * ri2 * scale
870 <             vpair = vpair + swi * vterm
871 <             epot = epot + vterm
870 >             vpair = vpair + vterm
871 >             epot = epot + sw*vterm
872              
873 <             dudx = dudx + pref * ri3 * ( uz_i(1) - 3.0d0 * ct_i * xhat*sc2)
874 <             dudy = dudy + pref * ri3 * ( uz_i(2) - 3.0d0 * ct_i * yhat*sc2)
875 <             dudz = dudz + pref * ri3 * ( uz_i(3) - 3.0d0 * ct_i * zhat*sc2)
873 >             dudx = dudx + sw*pref * ri3 * ( uz_i(1) - 3.0d0 * ct_i * xhat*sc2)
874 >             dudy = dudy + sw*pref * ri3 * ( uz_i(2) - 3.0d0 * ct_i * yhat*sc2)
875 >             dudz = dudz + sw*pref * ri3 * ( uz_i(3) - 3.0d0 * ct_i * zhat*sc2)
876              
877 <             duduz_i(1) = duduz_i(1) + pref * ri2 * xhat * scale
878 <             duduz_i(2) = duduz_i(2) + pref * ri2 * yhat * scale
879 <             duduz_i(3) = duduz_i(3) + pref * ri2 * zhat * scale
877 >             duduz_i(1) = duduz_i(1) + sw*pref * ri2 * xhat * scale
878 >             duduz_i(2) = duduz_i(2) + sw*pref * ri2 * yhat * scale
879 >             duduz_i(3) = duduz_i(3) + sw*pref * ri2 * zhat * scale
880            endif
881         endif
882 <
882 >      
883         if (j_is_Dipole) then
884  
885 <          pref = sw * pre22 * mu_i * mu_j
857 <
858 <          if (corrMethod .eq. 1) then
885 >          if (summationMethod .eq. UNDAMPED_WOLF) then
886               ri2 = riji * riji
887               ri3 = ri2 * riji
888               ri4 = ri2 * ri2
889  
890 +             pref = pre22 * mu_i * mu_j
891               vterm = pref * (ri3 - rcuti3) * (ct_ij - 3.0d0 * ct_i * ct_j)
892 <             vpair = vpair + swi * vterm
893 <             epot = epot + vterm
892 >             vpair = vpair + vterm
893 >             epot = epot + sw*vterm
894              
895               a1 = 5.0d0 * ct_i * ct_j - ct_ij
896              
897 <             dudx = dudx + pref*3.0d0*ri4 &
898 <                  *(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1)) - &
899 <                  pref*3.0d0*rcuti4*(a1*rcuti*d(1)-ct_i*uz_j(1)-ct_j*uz_i(1))
900 <             dudy = dudy + pref*3.0d0*ri4 &
901 <                  *(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2)) - &
902 <                  pref*3.0d0*rcuti4*(a1*rcuti*d(2)-ct_i*uz_j(2)-ct_j*uz_i(2))
903 <             dudz = dudz + pref*3.0d0*ri4 &
904 <                  *(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3)) - &
905 <                  pref*3.0d0*rcuti4*(a1*rcuti*d(3)-ct_i*uz_j(3)-ct_j*uz_i(3))
897 >             dudx = dudx + sw*pref*3.0d0*ri4 &
898 >                             * (a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1)) &
899 >                         - sw*pref*3.0d0*rcuti4 &
900 >                             * (a1*rcuti*d(1)-ct_i*uz_j(1)-ct_j*uz_i(1))
901 >             dudy = dudy + sw*pref*3.0d0*ri4 &
902 >                             * (a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2)) &
903 >                         - sw*pref*3.0d0*rcuti4 &
904 >                             * (a1*rcuti*d(2)-ct_i*uz_j(2)-ct_j*uz_i(2))
905 >             dudz = dudz + sw*pref*3.0d0*ri4 &
906 >                             * (a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3)) &
907 >                         - sw*pref*3.0d0*rcuti4 &
908 >                             * (a1*rcuti*d(3)-ct_i*uz_j(3)-ct_j*uz_i(3))
909              
910 <             duduz_i(1) = duduz_i(1) + pref*(ri3*(uz_j(1) - 3.0d0*ct_j*xhat) &
910 >             duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(uz_j(1)-3.0d0*ct_j*xhat) &
911                    - rcuti3*(uz_j(1) - 3.0d0*ct_j*d(1)*rcuti))
912 <             duduz_i(2) = duduz_i(2) + pref*(ri3*(uz_j(2) - 3.0d0*ct_j*yhat) &
912 >             duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(uz_j(2)-3.0d0*ct_j*yhat) &
913                    - rcuti3*(uz_j(2) - 3.0d0*ct_j*d(2)*rcuti))
914 <             duduz_i(3) = duduz_i(3) + pref*(ri3*(uz_j(3) - 3.0d0*ct_j*zhat) &
914 >             duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(uz_j(3)-3.0d0*ct_j*zhat) &
915                    - rcuti3*(uz_j(3) - 3.0d0*ct_j*d(3)*rcuti))
916 <             duduz_j(1) = duduz_j(1) + pref*(ri3*(uz_i(1) - 3.0d0*ct_i*xhat) &
916 >             duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(uz_i(1)-3.0d0*ct_i*xhat) &
917                    - rcuti3*(uz_i(1) - 3.0d0*ct_i*d(1)*rcuti))
918 <             duduz_j(2) = duduz_j(2) + pref*(ri3*(uz_i(2) - 3.0d0*ct_i*yhat) &
918 >             duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(uz_i(2)-3.0d0*ct_i*yhat) &
919                    - rcuti3*(uz_i(2) - 3.0d0*ct_i*d(2)*rcuti))
920 <             duduz_j(3) = duduz_j(3) + pref*(ri3*(uz_i(3) - 3.0d0*ct_i*zhat) &
920 >             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(uz_i(3)-3.0d0*ct_i*zhat) &
921                    - rcuti3*(uz_i(3) - 3.0d0*ct_i*d(3)*rcuti))
922 +
923            else
892 –            
924               if (i_is_SplitDipole) then
925                  if (j_is_SplitDipole) then
926                     BigR = sqrt(r2 + 0.25_dp * d_i * d_i + 0.25_dp * d_j * d_j)
# Line 916 | Line 947 | contains
947               ri4 = ri2 * ri2
948               sc2 = scale * scale
949              
950 +             pref = pre22 * mu_i * mu_j
951               vterm = pref * ri3 * (ct_ij - 3.0d0 * ct_i * ct_j * sc2)
952 <             vpair = vpair + swi * vterm
953 <             epot = epot + vterm
952 >             vpair = vpair + vterm
953 >             epot = epot + sw*vterm
954              
955               a1 = 5.0d0 * ct_i * ct_j * sc2 - ct_ij
956              
957 <             dudx = dudx + pref*3.0d0*ri4*scale &
958 <                  *(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
959 <             dudy = dudy + pref*3.0d0*ri4*scale &
960 <                  *(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
961 <             dudz = dudz + pref*3.0d0*ri4*scale &
962 <                  *(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
957 >             dudx = dudx + sw*pref*3.0d0*ri4*scale &
958 >                             *(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
959 >             dudy = dudy + sw*pref*3.0d0*ri4*scale &
960 >                             *(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
961 >             dudz = dudz + sw*pref*3.0d0*ri4*scale &
962 >                             *(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
963              
964 <             duduz_i(1) = duduz_i(1) + pref*ri3 &
965 <                  *(uz_j(1) - 3.0d0*ct_j*xhat*sc2)
966 <             duduz_i(2) = duduz_i(2) + pref*ri3 &
967 <                  *(uz_j(2) - 3.0d0*ct_j*yhat*sc2)
968 <             duduz_i(3) = duduz_i(3) + pref*ri3 &
969 <                  *(uz_j(3) - 3.0d0*ct_j*zhat*sc2)
964 >             duduz_i(1) = duduz_i(1) + sw*pref*ri3 &
965 >                                         *(uz_j(1) - 3.0d0*ct_j*xhat*sc2)
966 >             duduz_i(2) = duduz_i(2) + sw*pref*ri3 &
967 >                                         *(uz_j(2) - 3.0d0*ct_j*yhat*sc2)
968 >             duduz_i(3) = duduz_i(3) + sw*pref*ri3 &
969 >                                         *(uz_j(3) - 3.0d0*ct_j*zhat*sc2)
970              
971 <             duduz_j(1) = duduz_j(1) + pref*ri3 &
972 <                  *(uz_i(1) - 3.0d0*ct_i*xhat*sc2)
973 <             duduz_j(2) = duduz_j(2) + pref*ri3 &
974 <                  *(uz_i(2) - 3.0d0*ct_i*yhat*sc2)
975 <             duduz_j(3) = duduz_j(3) + pref*ri3 &
976 <                  *(uz_i(3) - 3.0d0*ct_i*zhat*sc2)
971 >             duduz_j(1) = duduz_j(1) + sw*pref*ri3 &
972 >                                         *(uz_i(1) - 3.0d0*ct_i*xhat*sc2)
973 >             duduz_j(2) = duduz_j(2) + sw*pref*ri3 &
974 >                                         *(uz_i(2) - 3.0d0*ct_i*yhat*sc2)
975 >             duduz_j(3) = duduz_j(3) + sw*pref*ri3 &
976 >                                         *(uz_i(3) - 3.0d0*ct_i*zhat*sc2)
977            endif
978         endif
979      endif
# Line 956 | Line 988 | contains
988            cy2 = cy_i * cy_i
989            cz2 = cz_i * cz_i
990  
991 <          pref = sw * pre14 * q_j / 3.0_dp
992 <
961 <          if (corrMethod .eq. 1) then
991 >          if (summationMethod .eq. UNDAMPED_WOLF) then
992 >             pref = pre14 * q_j / 3.0_dp
993               vterm1 = pref * ri3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
994                    qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
995                    qzz_i * (3.0_dp*cz2 - 1.0_dp) )
996               vterm2 = pref * rcuti3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
997                    qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
998                    qzz_i * (3.0_dp*cz2 - 1.0_dp) )
999 <             vpair = vpair + swi * ( vterm1 - vterm2 )
1000 <             epot = epot + ( vterm1 - vterm2 )
999 >             vpair = vpair + ( vterm1 - vterm2 )
1000 >             epot = epot + sw*( vterm1 - vterm2 )
1001              
1002 <             dudx = dudx - (5.0_dp*(vterm1*riji*xhat - vterm2*rcuti2*d(1))) + &
1003 <                  pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(1)) - &
1002 >             dudx = dudx - sw*(5.0_dp*(vterm1*riji*xhat-vterm2*rcuti2*d(1))) +&
1003 >                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(1)) - &
1004                    qxx_i*2.0_dp*(xhat - rcuti*d(1))) + &
1005                    (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(1)) - &
1006                    qyy_i*2.0_dp*(xhat - rcuti*d(1))) + &
1007                    (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(1)) - &
1008                    qzz_i*2.0_dp*(xhat - rcuti*d(1))) )
1009 <             dudy = dudy - (5.0_dp*(vterm1*riji*yhat - vterm2*rcuti2*d(2))) + &
1010 <                  pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(2)) - &
1009 >             dudy = dudy - sw*(5.0_dp*(vterm1*riji*yhat-vterm2*rcuti2*d(2))) +&
1010 >                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(2)) - &
1011                    qxx_i*2.0_dp*(yhat - rcuti*d(2))) + &
1012                    (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(2)) - &
1013                    qyy_i*2.0_dp*(yhat - rcuti*d(2))) + &
1014                    (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(2)) - &
1015                    qzz_i*2.0_dp*(yhat - rcuti*d(2))) )
1016 <             dudz = dudz - (5.0_dp*(vterm1*riji*zhat - vterm2*rcuti2*d(3))) + &
1017 <                  pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(3)) - &
1016 >             dudz = dudz - sw*(5.0_dp*(vterm1*riji*zhat-vterm2*rcuti2*d(3))) +&
1017 >                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(3)) - &
1018                    qxx_i*2.0_dp*(zhat - rcuti*d(3))) + &
1019                    (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(3)) - &
1020                    qyy_i*2.0_dp*(zhat - rcuti*d(3))) + &
1021                    (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(3)) - &
1022                    qzz_i*2.0_dp*(zhat - rcuti*d(3))) )
1023              
1024 <             dudux_i(1) = dudux_i(1) + pref * (ri3*(qxx_i*6.0_dp*cx_i*xhat) - &
1024 >             dudux_i(1) = dudux_i(1) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*xhat) -&
1025                    rcuti4*(qxx_i*6.0_dp*cx_i*d(1)))
1026 <             dudux_i(2) = dudux_i(2) + pref * (ri3*(qxx_i*6.0_dp*cx_i*yhat) - &
1026 >             dudux_i(2) = dudux_i(2) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*yhat) -&
1027                    rcuti4*(qxx_i*6.0_dp*cx_i*d(2)))
1028 <             dudux_i(3) = dudux_i(3) + pref * (ri3*(qxx_i*6.0_dp*cx_i*zhat) - &
1028 >             dudux_i(3) = dudux_i(3) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*zhat) -&
1029                    rcuti4*(qxx_i*6.0_dp*cx_i*d(3)))
1030              
1031 <             duduy_i(1) = duduy_i(1) + pref * (ri3*(qyy_i*6.0_dp*cy_i*xhat) - &
1031 >             duduy_i(1) = duduy_i(1) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*xhat) -&
1032                    rcuti4*(qyy_i*6.0_dp*cx_i*d(1)))
1033 <             duduy_i(2) = duduy_i(2) + pref * (ri3*(qyy_i*6.0_dp*cy_i*yhat) - &
1033 >             duduy_i(2) = duduy_i(2) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*yhat) -&
1034                    rcuti4*(qyy_i*6.0_dp*cx_i*d(2)))
1035 <             duduy_i(3) = duduy_i(3) + pref * (ri3*(qyy_i*6.0_dp*cy_i*zhat) - &
1035 >             duduy_i(3) = duduy_i(3) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*zhat) -&
1036                    rcuti4*(qyy_i*6.0_dp*cx_i*d(3)))
1037              
1038 <             duduz_i(1) = duduz_i(1) + pref * (ri3*(qzz_i*6.0_dp*cz_i*xhat) - &
1038 >             duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*xhat) -&
1039                    rcuti4*(qzz_i*6.0_dp*cx_i*d(1)))
1040 <             duduz_i(2) = duduz_i(2) + pref * (ri3*(qzz_i*6.0_dp*cz_i*yhat) - &
1040 >             duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*yhat) -&
1041                    rcuti4*(qzz_i*6.0_dp*cx_i*d(2)))
1042 <             duduz_i(3) = duduz_i(3) + pref * (ri3*(qzz_i*6.0_dp*cz_i*zhat) - &
1042 >             duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*zhat) -&
1043                    rcuti4*(qzz_i*6.0_dp*cx_i*d(3)))
1044  
1045            else
1046 +             pref = pre14 * q_j / 3.0_dp
1047               vterm = pref * ri3 * (qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1048                    qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1049                    qzz_i * (3.0_dp*cz2 - 1.0_dp))
1050 <             vpair = vpair + swi * vterm
1051 <             epot = epot + vterm
1050 >             vpair = vpair + vterm
1051 >             epot = epot + sw*vterm
1052              
1053 <             dudx = dudx - 5.0_dp*vterm*riji*xhat + pref * ri4 * ( &
1053 >             dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref*ri4 * ( &
1054                    qxx_i*(6.0_dp*cx_i*ux_i(1) - 2.0_dp*xhat) + &
1055                    qyy_i*(6.0_dp*cy_i*uy_i(1) - 2.0_dp*xhat) + &
1056                    qzz_i*(6.0_dp*cz_i*uz_i(1) - 2.0_dp*xhat) )
1057 <             dudy = dudy - 5.0_dp*vterm*riji*yhat + pref * ri4 * ( &
1057 >             dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref*ri4 * ( &
1058                    qxx_i*(6.0_dp*cx_i*ux_i(2) - 2.0_dp*yhat) + &
1059                    qyy_i*(6.0_dp*cy_i*uy_i(2) - 2.0_dp*yhat) + &
1060                    qzz_i*(6.0_dp*cz_i*uz_i(2) - 2.0_dp*yhat) )
1061 <             dudz = dudz - 5.0_dp*vterm*riji*zhat + pref * ri4 * ( &
1061 >             dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref*ri4 * ( &
1062                    qxx_i*(6.0_dp*cx_i*ux_i(3) - 2.0_dp*zhat) + &
1063                    qyy_i*(6.0_dp*cy_i*uy_i(3) - 2.0_dp*zhat) + &
1064                    qzz_i*(6.0_dp*cz_i*uz_i(3) - 2.0_dp*zhat) )
1065              
1066 <             dudux_i(1) = dudux_i(1) + pref * ri3*(qxx_i*6.0_dp*cx_i*xhat)
1067 <             dudux_i(2) = dudux_i(2) + pref * ri3*(qxx_i*6.0_dp*cx_i*yhat)
1068 <             dudux_i(3) = dudux_i(3) + pref * ri3*(qxx_i*6.0_dp*cx_i*zhat)
1066 >             dudux_i(1) = dudux_i(1) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*xhat)
1067 >             dudux_i(2) = dudux_i(2) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*yhat)
1068 >             dudux_i(3) = dudux_i(3) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*zhat)
1069              
1070 <             duduy_i(1) = duduy_i(1) + pref * ri3*(qyy_i*6.0_dp*cy_i*xhat)
1071 <             duduy_i(2) = duduy_i(2) + pref * ri3*(qyy_i*6.0_dp*cy_i*yhat)
1072 <             duduy_i(3) = duduy_i(3) + pref * ri3*(qyy_i*6.0_dp*cy_i*zhat)
1070 >             duduy_i(1) = duduy_i(1) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*xhat)
1071 >             duduy_i(2) = duduy_i(2) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*yhat)
1072 >             duduy_i(3) = duduy_i(3) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*zhat)
1073              
1074 <             duduz_i(1) = duduz_i(1) + pref * ri3*(qzz_i*6.0_dp*cz_i*xhat)
1075 <             duduz_i(2) = duduz_i(2) + pref * ri3*(qzz_i*6.0_dp*cz_i*yhat)
1076 <             duduz_i(3) = duduz_i(3) + pref * ri3*(qzz_i*6.0_dp*cz_i*zhat)
1074 >             duduz_i(1) = duduz_i(1) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*xhat)
1075 >             duduz_i(2) = duduz_i(2) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*yhat)
1076 >             duduz_i(3) = duduz_i(3) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*zhat)
1077            endif
1078         endif
1079      endif

Diff Legend

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