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 2331 by chuckv, Mon Sep 26 18:38:02 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 +  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 <
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 +    if (summationMethod .eq. UNDAMPED_WOLF) is_Undamped_Wolf = .true.
140 +    if (summationMethod .eq. DAMPED_WOLF) is_Damped_Wolf = .true.
141    end subroutine setElectrostaticSummationMethod
142  
143    subroutine setElectrostaticCutoffRadius(thisRcut)
# Line 357 | Line 370 | contains
370  
371    subroutine checkSummationMethod()
372  
373 +    if (.not.haveDefaultCutoff) then
374 +       call handleError("checkSummationMethod", "no Default Cutoff set!")
375 +    endif
376 +
377 +    rcuti = 1.0d0 / defaultCutoff
378 +    rcuti2 = rcuti*rcuti
379 +    rcuti3 = rcuti2*rcuti
380 +    rcuti4 = rcuti2*rcuti2
381 +
382      if (summationMethod .eq. DAMPED_WOLF) then
383         if (.not.haveDWAconstants) then
384            
# Line 369 | Line 391 | contains
391            endif
392  
393            constEXP = exp(-dampingAlpha*dampingAlpha*defaultCutoff*defaultCutoff)
394 <          constERFC = erfc(dampingAlpha*defaultCutoff)
394 >          constERFC = derfc(dampingAlpha*defaultCutoff)
395            
396            haveDWAconstants = .true.
397         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
429 <    real (kind=dp) :: rcuti2, rcuti3, rcuti4
449 >    real (kind=dp) :: scale, sc2, bigR
450  
451      if (.not.allocated(ElectrostaticMap)) then
452         call handleError("electrostatic", "no ElectrostaticMap was present before first call of do_electrostatic_pair!")
# Line 435 | Line 455 | contains
455  
456      if (.not.summationMethodChecked) then
457         call checkSummationMethod()
458 +      
459      endif
460  
461  
# Line 454 | Line 475 | contains
475      yhat = d(2) * riji
476      zhat = d(3) * riji
477  
457 –    rcuti2 = rcuti*rcuti
458 –    rcuti3 = rcuti2*rcuti
459 –    rcuti4 = rcuti2*rcuti2
460 –
461 –    swi = 1.0d0 / sw
462 –
478      !! logicals
479      i_is_Charge = ElectrostaticMap(me1)%is_Charge
480      i_is_Dipole = ElectrostaticMap(me1)%is_Dipole
# Line 578 | Line 593 | contains
593         cz_j = uz_j(1)*xhat + uz_j(2)*yhat + uz_j(3)*zhat
594      endif
595    
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 –
596      epot = 0.0_dp
597      dudx = 0.0_dp
598      dudy = 0.0_dp
# Line 603 | Line 610 | contains
610  
611         if (j_is_Charge) then
612  
613 <          if (corrMethod .eq. 1) then
607 <             vterm = pre11 * q_i * q_j * (riji - rcuti)
613 >          if (summationMethod .eq. UNDAMPED_WOLF) then
614  
615 +             vterm = pre11 * q_i * q_j * (riji - rcuti)
616               vpair = vpair + vterm
617 <             epot = epot + sw * vterm
617 >             epot = epot + sw*vterm
618              
619               dudr  = - sw * pre11 * q_i * q_j * (riji*riji*riji - rcuti2*rcuti)
620              
# Line 616 | Line 623 | contains
623               dudz = dudz + dudr * d(3)
624  
625            else
619 –             vterm = pre11 * q_i * q_j * riji
626  
627 +             vterm = pre11 * q_i * q_j * riji
628               vpair = vpair + vterm
629 <             epot = epot + sw * vterm
629 >             epot = epot + sw*vterm
630              
631               dudr  = - sw * vterm * riji
632              
# Line 633 | Line 640 | contains
640  
641         if (j_is_Dipole) then
642  
643 <          pref = sw * pre12 * q_i * mu_j
643 >          pref = pre12 * q_i * mu_j
644  
645 <          if (corrMethod .eq. 1) then
645 >          if (summationMethod .eq. UNDAMPED_WOLF) then
646               ri2 = riji * riji
647               ri3 = ri2 * riji
648  
649 +             pref = pre12 * q_i * mu_j
650               vterm = - pref * ct_j * (ri2 - rcuti2)
651 <             vpair = vpair + swi*vterm
652 <             epot = epot + vterm
651 >             vpair = vpair + vterm
652 >             epot = epot + sw*vterm
653              
654               !! this has a + sign in the () because the rij vector is
655               !! r_j - r_i and the charge-dipole potential takes the origin
656               !! as the point dipole, which is atom j in this case.
657              
658 <             dudx = dudx - pref * ( ri3*( uz_j(1) - 3.0d0*ct_j*xhat) &
658 >             dudx = dudx - sw*pref * ( ri3*( uz_j(1) - 3.0d0*ct_j*xhat) &
659                    - rcuti3*( uz_j(1) - 3.0d0*ct_j*d(1)*rcuti ) )
660 <             dudy = dudy - pref * ( ri3*( uz_j(2) - 3.0d0*ct_j*yhat) &
660 >             dudy = dudy - sw*pref * ( ri3*( uz_j(2) - 3.0d0*ct_j*yhat) &
661                    - rcuti3*( uz_j(2) - 3.0d0*ct_j*d(2)*rcuti ) )
662 <             dudz = dudz - pref * ( ri3*( uz_j(3) - 3.0d0*ct_j*zhat) &
662 >             dudz = dudz - sw*pref * ( ri3*( uz_j(3) - 3.0d0*ct_j*zhat) &
663                    - rcuti3*( uz_j(3) - 3.0d0*ct_j*d(3)*rcuti ) )
664              
665 <             duduz_j(1) = duduz_j(1) - pref*( ri2*xhat - d(1)*rcuti3 )
666 <             duduz_j(2) = duduz_j(2) - pref*( ri2*yhat - d(2)*rcuti3 )
667 <             duduz_j(3) = duduz_j(3) - pref*( ri2*zhat - d(3)*rcuti3 )
665 >             duduz_j(1) = duduz_j(1) - sw*pref*( ri2*xhat - d(1)*rcuti3 )
666 >             duduz_j(2) = duduz_j(2) - sw*pref*( ri2*yhat - d(2)*rcuti3 )
667 >             duduz_j(3) = duduz_j(3) - sw*pref*( ri2*zhat - d(3)*rcuti3 )
668  
669            else
670               if (j_is_SplitDipole) then
# Line 671 | Line 679 | contains
679               ri2 = ri * ri
680               ri3 = ri2 * ri
681               sc2 = scale * scale
682 <            
682 >
683 >             pref = pre12 * q_i * mu_j
684               vterm = - pref * ct_j * ri2 * scale
685 <             vpair = vpair + swi * vterm
686 <             epot = epot + vterm
685 >             vpair = vpair + vterm
686 >             epot = epot + sw*vterm
687              
688               !! this has a + sign in the () because the rij vector is
689               !! r_j - r_i and the charge-dipole potential takes the origin
690               !! as the point dipole, which is atom j in this case.
691              
692 <             dudx = dudx - pref * ri3 * ( uz_j(1) - 3.0d0*ct_j*xhat*sc2)
693 <             dudy = dudy - pref * ri3 * ( uz_j(2) - 3.0d0*ct_j*yhat*sc2)
694 <             dudz = dudz - pref * ri3 * ( uz_j(3) - 3.0d0*ct_j*zhat*sc2)
692 >             dudx = dudx - sw*pref * ri3 * ( uz_j(1) - 3.0d0*ct_j*xhat*sc2)
693 >             dudy = dudy - sw*pref * ri3 * ( uz_j(2) - 3.0d0*ct_j*yhat*sc2)
694 >             dudz = dudz - sw*pref * ri3 * ( uz_j(3) - 3.0d0*ct_j*zhat*sc2)
695              
696 <             duduz_j(1) = duduz_j(1) - pref * ri2 * xhat * scale
697 <             duduz_j(2) = duduz_j(2) - pref * ri2 * yhat * scale
698 <             duduz_j(3) = duduz_j(3) - pref * ri2 * zhat * scale
696 >             duduz_j(1) = duduz_j(1) - sw*pref * ri2 * xhat * scale
697 >             duduz_j(2) = duduz_j(2) - sw*pref * ri2 * yhat * scale
698 >             duduz_j(3) = duduz_j(3) - sw*pref * ri2 * zhat * scale
699  
700            endif
701         endif
# Line 699 | Line 708 | contains
708            cy2 = cy_j * cy_j
709            cz2 = cz_j * cz_j
710  
711 <
712 <          pref =  sw * pre14 * q_i / 3.0_dp
704 <
705 <          if (corrMethod .eq. 1) then
711 >          if (summationMethod .eq. UNDAMPED_WOLF) then
712 >             pref =  pre14 * q_i / 3.0_dp
713               vterm1 = pref * ri3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
714                    qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
715                    qzz_j * (3.0_dp*cz2 - 1.0_dp) )
716               vterm2 = pref * rcuti3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
717                    qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
718                    qzz_j * (3.0_dp*cz2 - 1.0_dp) )
719 <             vpair = vpair + swi*( vterm1 - vterm2 )
720 <             epot = epot + ( vterm1 - vterm2 )
719 >             vpair = vpair + ( vterm1 - vterm2 )
720 >             epot = epot + sw*( vterm1 - vterm2 )
721              
722               dudx = dudx - (5.0_dp * &
723 <                  (vterm1*riji*xhat - vterm2*rcuti2*d(1))) + pref * ( &
723 >                  (vterm1*riji*xhat - vterm2*rcuti2*d(1))) + sw*pref * ( &
724                    (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(1)) - &
725                    qxx_j*2.0_dp*(xhat - rcuti*d(1))) + &
726                    (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(1)) - &
# Line 721 | Line 728 | contains
728                    (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(1)) - &
729                    qzz_j*2.0_dp*(xhat - rcuti*d(1))) )
730               dudy = dudy - (5.0_dp * &
731 <                  (vterm1*riji*yhat - vterm2*rcuti2*d(2))) + pref * ( &
731 >                  (vterm1*riji*yhat - vterm2*rcuti2*d(2))) + sw*pref * ( &
732                    (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(2)) - &
733                    qxx_j*2.0_dp*(yhat - rcuti*d(2))) + &
734                    (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(2)) - &
# Line 729 | Line 736 | contains
736                    (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(2)) - &
737                    qzz_j*2.0_dp*(yhat - rcuti*d(2))) )
738               dudz = dudz - (5.0_dp * &
739 <                  (vterm1*riji*zhat - vterm2*rcuti2*d(3))) + pref * ( &
739 >                  (vterm1*riji*zhat - vterm2*rcuti2*d(3))) + sw*pref * ( &
740                    (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(3)) - &
741                    qxx_j*2.0_dp*(zhat - rcuti*d(3))) + &
742                    (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(3)) - &
# Line 737 | Line 744 | contains
744                    (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(3)) - &
745                    qzz_j*2.0_dp*(zhat - rcuti*d(3))) )
746              
747 <             dudux_j(1) = dudux_j(1) + pref * (ri3*(qxx_j*6.0_dp*cx_j*xhat) - &
747 >             dudux_j(1) = dudux_j(1) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*xhat) -&
748                    rcuti4*(qxx_j*6.0_dp*cx_j*d(1)))
749 <             dudux_j(2) = dudux_j(2) + pref * (ri3*(qxx_j*6.0_dp*cx_j*yhat) - &
749 >             dudux_j(2) = dudux_j(2) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*yhat) -&
750                    rcuti4*(qxx_j*6.0_dp*cx_j*d(2)))
751 <             dudux_j(3) = dudux_j(3) + pref * (ri3*(qxx_j*6.0_dp*cx_j*zhat) - &
751 >             dudux_j(3) = dudux_j(3) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*zhat) -&
752                    rcuti4*(qxx_j*6.0_dp*cx_j*d(3)))
753              
754 <             duduy_j(1) = duduy_j(1) + pref * (ri3*(qyy_j*6.0_dp*cy_j*xhat) - &
754 >             duduy_j(1) = duduy_j(1) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*xhat) -&
755                    rcuti4*(qyy_j*6.0_dp*cx_j*d(1)))
756 <             duduy_j(2) = duduy_j(2) + pref * (ri3*(qyy_j*6.0_dp*cy_j*yhat) - &
756 >             duduy_j(2) = duduy_j(2) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*yhat) -&
757                    rcuti4*(qyy_j*6.0_dp*cx_j*d(2)))
758 <             duduy_j(3) = duduy_j(3) + pref * (ri3*(qyy_j*6.0_dp*cy_j*zhat) - &
758 >             duduy_j(3) = duduy_j(3) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*zhat) -&
759                    rcuti4*(qyy_j*6.0_dp*cx_j*d(3)))
760              
761 <             duduz_j(1) = duduz_j(1) + pref * (ri3*(qzz_j*6.0_dp*cz_j*xhat) - &
761 >             duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*xhat) -&
762                    rcuti4*(qzz_j*6.0_dp*cx_j*d(1)))
763 <             duduz_j(2) = duduz_j(2) + pref * (ri3*(qzz_j*6.0_dp*cz_j*yhat) - &
763 >             duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*yhat) -&
764                    rcuti4*(qzz_j*6.0_dp*cx_j*d(2)))
765 <             duduz_j(3) = duduz_j(3) + pref * (ri3*(qzz_j*6.0_dp*cz_j*zhat) - &
765 >             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*zhat) -&
766                    rcuti4*(qzz_j*6.0_dp*cx_j*d(3)))
767          
768            else
769 +             pref =  pre14 * q_i / 3.0_dp
770               vterm = pref * ri3 * (qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
771                    qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
772                    qzz_j * (3.0_dp*cz2 - 1.0_dp))
773 <             vpair = vpair + swi * vterm
774 <             epot = epot + vterm
773 >             vpair = vpair + vterm
774 >             epot = epot + sw*vterm
775              
776 <             dudx = dudx - 5.0_dp*vterm*riji*xhat + pref * ri4 * ( &
776 >             dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref * ri4 * ( &
777                    qxx_j*(6.0_dp*cx_j*ux_j(1) - 2.0_dp*xhat) + &
778                    qyy_j*(6.0_dp*cy_j*uy_j(1) - 2.0_dp*xhat) + &
779                    qzz_j*(6.0_dp*cz_j*uz_j(1) - 2.0_dp*xhat) )
780 <             dudy = dudy - 5.0_dp*vterm*riji*yhat + pref * ri4 * ( &
780 >             dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref * ri4 * ( &
781                    qxx_j*(6.0_dp*cx_j*ux_j(2) - 2.0_dp*yhat) + &
782                    qyy_j*(6.0_dp*cy_j*uy_j(2) - 2.0_dp*yhat) + &
783                    qzz_j*(6.0_dp*cz_j*uz_j(2) - 2.0_dp*yhat) )
784 <             dudz = dudz - 5.0_dp*vterm*riji*zhat + pref * ri4 * ( &
784 >             dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref * ri4 * ( &
785                    qxx_j*(6.0_dp*cx_j*ux_j(3) - 2.0_dp*zhat) + &
786                    qyy_j*(6.0_dp*cy_j*uy_j(3) - 2.0_dp*zhat) + &
787                    qzz_j*(6.0_dp*cz_j*uz_j(3) - 2.0_dp*zhat) )
788              
789 <             dudux_j(1) = dudux_j(1) + pref * ri3*(qxx_j*6.0_dp*cx_j*xhat)
790 <             dudux_j(2) = dudux_j(2) + pref * ri3*(qxx_j*6.0_dp*cx_j*yhat)
791 <             dudux_j(3) = dudux_j(3) + pref * ri3*(qxx_j*6.0_dp*cx_j*zhat)
789 >             dudux_j(1) = dudux_j(1) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*xhat)
790 >             dudux_j(2) = dudux_j(2) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*yhat)
791 >             dudux_j(3) = dudux_j(3) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*zhat)
792              
793 <             duduy_j(1) = duduy_j(1) + pref * ri3*(qyy_j*6.0_dp*cy_j*xhat)
794 <             duduy_j(2) = duduy_j(2) + pref * ri3*(qyy_j*6.0_dp*cy_j*yhat)
795 <             duduy_j(3) = duduy_j(3) + pref * ri3*(qyy_j*6.0_dp*cy_j*zhat)
793 >             duduy_j(1) = duduy_j(1) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*xhat)
794 >             duduy_j(2) = duduy_j(2) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*yhat)
795 >             duduy_j(3) = duduy_j(3) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*zhat)
796              
797 <             duduz_j(1) = duduz_j(1) + pref * ri3*(qzz_j*6.0_dp*cz_j*xhat)
798 <             duduz_j(2) = duduz_j(2) + pref * ri3*(qzz_j*6.0_dp*cz_j*yhat)
799 <             duduz_j(3) = duduz_j(3) + pref * ri3*(qzz_j*6.0_dp*cz_j*zhat)
797 >             duduz_j(1) = duduz_j(1) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*xhat)
798 >             duduz_j(2) = duduz_j(2) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*yhat)
799 >             duduz_j(3) = duduz_j(3) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*zhat)
800            
801            endif
802         endif
# Line 797 | Line 805 | contains
805      if (i_is_Dipole) then
806  
807         if (j_is_Charge) then
808 <
809 <          pref = sw * pre12 * q_j * mu_i
810 <
811 <          if (corrMethod .eq. 1) then
808 >          
809 >          pref = pre12 * q_j * mu_i
810 >          
811 >          if (summationMethod .eq. UNDAMPED_WOLF) then
812               ri2 = riji * riji
813               ri3 = ri2 * riji
814  
815 +             pref = pre12 * q_j * mu_i
816               vterm = pref * ct_i * (ri2 - rcuti2)
817 <             vpair = vpair + swi * vterm
818 <             epot = epot + vterm
817 >             vpair = vpair + vterm
818 >             epot = epot + sw*vterm
819              
820               !! this has a + sign in the () because the rij vector is
821               !! r_j - r_i and the charge-dipole potential takes the origin
822               !! as the point dipole, which is atom j in this case.
823              
824 <             dudx = dudx + pref * ( ri3*( uz_i(1) - 3.0d0*ct_i*xhat) &
824 >             dudx = dudx + sw*pref * ( ri3*( uz_i(1) - 3.0d0*ct_i*xhat) &
825                    - rcuti3*( uz_i(1) - 3.0d0*ct_i*d(1)*rcuti ) )
826 <             dudy = dudy + pref * ( ri3*( uz_i(2) - 3.0d0*ct_i*yhat) &
826 >             dudy = dudy + sw*pref * ( ri3*( uz_i(2) - 3.0d0*ct_i*yhat) &
827                    - rcuti3*( uz_i(2) - 3.0d0*ct_i*d(2)*rcuti ) )
828 <             dudz = dudz + pref * ( ri3*( uz_i(3) - 3.0d0*ct_i*zhat) &
828 >             dudz = dudz + sw*pref * ( ri3*( uz_i(3) - 3.0d0*ct_i*zhat) &
829                    - rcuti3*( uz_i(3) - 3.0d0*ct_i*d(3)*rcuti ) )
830              
831 <             duduz_i(1) = duduz_i(1) - pref*( ri2*xhat - d(1)*rcuti3 )
832 <             duduz_i(2) = duduz_i(2) - pref*( ri2*yhat - d(2)*rcuti3 )
833 <             duduz_i(3) = duduz_i(3) - pref*( ri2*zhat - d(3)*rcuti3 )
831 >             duduz_i(1) = duduz_i(1) - sw*pref*( ri2*xhat - d(1)*rcuti3 )
832 >             duduz_i(2) = duduz_i(2) - sw*pref*( ri2*yhat - d(2)*rcuti3 )
833 >             duduz_i(3) = duduz_i(3) - sw*pref*( ri2*zhat - d(3)*rcuti3 )
834  
835            else
836               if (i_is_SplitDipole) then
# Line 836 | Line 845 | contains
845               ri2 = ri * ri
846               ri3 = ri2 * ri
847               sc2 = scale * scale
848 <            
848 >
849 >             pref = pre12 * q_j * mu_i
850               vterm = pref * ct_i * ri2 * scale
851 <             vpair = vpair + swi * vterm
852 <             epot = epot + vterm
851 >             vpair = vpair + vterm
852 >             epot = epot + sw*vterm
853              
854 <             dudx = dudx + pref * ri3 * ( uz_i(1) - 3.0d0 * ct_i * xhat*sc2)
855 <             dudy = dudy + pref * ri3 * ( uz_i(2) - 3.0d0 * ct_i * yhat*sc2)
856 <             dudz = dudz + pref * ri3 * ( uz_i(3) - 3.0d0 * ct_i * zhat*sc2)
854 >             dudx = dudx + sw*pref * ri3 * ( uz_i(1) - 3.0d0 * ct_i * xhat*sc2)
855 >             dudy = dudy + sw*pref * ri3 * ( uz_i(2) - 3.0d0 * ct_i * yhat*sc2)
856 >             dudz = dudz + sw*pref * ri3 * ( uz_i(3) - 3.0d0 * ct_i * zhat*sc2)
857              
858 <             duduz_i(1) = duduz_i(1) + pref * ri2 * xhat * scale
859 <             duduz_i(2) = duduz_i(2) + pref * ri2 * yhat * scale
860 <             duduz_i(3) = duduz_i(3) + pref * ri2 * zhat * scale
858 >             duduz_i(1) = duduz_i(1) + sw*pref * ri2 * xhat * scale
859 >             duduz_i(2) = duduz_i(2) + sw*pref * ri2 * yhat * scale
860 >             duduz_i(3) = duduz_i(3) + sw*pref * ri2 * zhat * scale
861            endif
862         endif
863 <
863 >      
864         if (j_is_Dipole) then
865  
866 <          pref = sw * pre22 * mu_i * mu_j
857 <
858 <          if (corrMethod .eq. 1) then
866 >          if (summationMethod .eq. UNDAMPED_WOLF) then
867               ri2 = riji * riji
868               ri3 = ri2 * riji
869               ri4 = ri2 * ri2
870  
871 +             pref = pre22 * mu_i * mu_j
872               vterm = pref * (ri3 - rcuti3) * (ct_ij - 3.0d0 * ct_i * ct_j)
873 <             vpair = vpair + swi * vterm
874 <             epot = epot + vterm
873 >             vpair = vpair + vterm
874 >             epot = epot + sw*vterm
875              
876               a1 = 5.0d0 * ct_i * ct_j - ct_ij
877              
878 <             dudx = dudx + pref*3.0d0*ri4 &
879 <                  *(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1)) - &
880 <                  pref*3.0d0*rcuti4*(a1*rcuti*d(1)-ct_i*uz_j(1)-ct_j*uz_i(1))
881 <             dudy = dudy + pref*3.0d0*ri4 &
882 <                  *(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2)) - &
883 <                  pref*3.0d0*rcuti4*(a1*rcuti*d(2)-ct_i*uz_j(2)-ct_j*uz_i(2))
884 <             dudz = dudz + pref*3.0d0*ri4 &
885 <                  *(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3)) - &
886 <                  pref*3.0d0*rcuti4*(a1*rcuti*d(3)-ct_i*uz_j(3)-ct_j*uz_i(3))
878 >             dudx = dudx + sw*pref*3.0d0*ri4 &
879 >                             * (a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1)) &
880 >                         - sw*pref*3.0d0*rcuti4 &
881 >                             * (a1*rcuti*d(1)-ct_i*uz_j(1)-ct_j*uz_i(1))
882 >             dudy = dudy + sw*pref*3.0d0*ri4 &
883 >                             * (a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2)) &
884 >                         - sw*pref*3.0d0*rcuti4 &
885 >                             * (a1*rcuti*d(2)-ct_i*uz_j(2)-ct_j*uz_i(2))
886 >             dudz = dudz + sw*pref*3.0d0*ri4 &
887 >                             * (a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3)) &
888 >                         - sw*pref*3.0d0*rcuti4 &
889 >                             * (a1*rcuti*d(3)-ct_i*uz_j(3)-ct_j*uz_i(3))
890              
891 <             duduz_i(1) = duduz_i(1) + pref*(ri3*(uz_j(1) - 3.0d0*ct_j*xhat) &
891 >             duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(uz_j(1)-3.0d0*ct_j*xhat) &
892                    - rcuti3*(uz_j(1) - 3.0d0*ct_j*d(1)*rcuti))
893 <             duduz_i(2) = duduz_i(2) + pref*(ri3*(uz_j(2) - 3.0d0*ct_j*yhat) &
893 >             duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(uz_j(2)-3.0d0*ct_j*yhat) &
894                    - rcuti3*(uz_j(2) - 3.0d0*ct_j*d(2)*rcuti))
895 <             duduz_i(3) = duduz_i(3) + pref*(ri3*(uz_j(3) - 3.0d0*ct_j*zhat) &
895 >             duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(uz_j(3)-3.0d0*ct_j*zhat) &
896                    - rcuti3*(uz_j(3) - 3.0d0*ct_j*d(3)*rcuti))
897 <             duduz_j(1) = duduz_j(1) + pref*(ri3*(uz_i(1) - 3.0d0*ct_i*xhat) &
897 >             duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(uz_i(1)-3.0d0*ct_i*xhat) &
898                    - rcuti3*(uz_i(1) - 3.0d0*ct_i*d(1)*rcuti))
899 <             duduz_j(2) = duduz_j(2) + pref*(ri3*(uz_i(2) - 3.0d0*ct_i*yhat) &
899 >             duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(uz_i(2)-3.0d0*ct_i*yhat) &
900                    - rcuti3*(uz_i(2) - 3.0d0*ct_i*d(2)*rcuti))
901 <             duduz_j(3) = duduz_j(3) + pref*(ri3*(uz_i(3) - 3.0d0*ct_i*zhat) &
901 >             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(uz_i(3)-3.0d0*ct_i*zhat) &
902                    - rcuti3*(uz_i(3) - 3.0d0*ct_i*d(3)*rcuti))
903 +
904            else
892 –            
905               if (i_is_SplitDipole) then
906                  if (j_is_SplitDipole) then
907                     BigR = sqrt(r2 + 0.25_dp * d_i * d_i + 0.25_dp * d_j * d_j)
# Line 916 | Line 928 | contains
928               ri4 = ri2 * ri2
929               sc2 = scale * scale
930              
931 +             pref = pre22 * mu_i * mu_j
932               vterm = pref * ri3 * (ct_ij - 3.0d0 * ct_i * ct_j * sc2)
933 <             vpair = vpair + swi * vterm
934 <             epot = epot + vterm
933 >             vpair = vpair + vterm
934 >             epot = epot + sw*vterm
935              
936               a1 = 5.0d0 * ct_i * ct_j * sc2 - ct_ij
937              
938 <             dudx = dudx + pref*3.0d0*ri4*scale &
939 <                  *(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
940 <             dudy = dudy + pref*3.0d0*ri4*scale &
941 <                  *(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
942 <             dudz = dudz + pref*3.0d0*ri4*scale &
943 <                  *(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
938 >             dudx = dudx + sw*pref*3.0d0*ri4*scale &
939 >                             *(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
940 >             dudy = dudy + sw*pref*3.0d0*ri4*scale &
941 >                             *(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
942 >             dudz = dudz + sw*pref*3.0d0*ri4*scale &
943 >                             *(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
944              
945 <             duduz_i(1) = duduz_i(1) + pref*ri3 &
946 <                  *(uz_j(1) - 3.0d0*ct_j*xhat*sc2)
947 <             duduz_i(2) = duduz_i(2) + pref*ri3 &
948 <                  *(uz_j(2) - 3.0d0*ct_j*yhat*sc2)
949 <             duduz_i(3) = duduz_i(3) + pref*ri3 &
950 <                  *(uz_j(3) - 3.0d0*ct_j*zhat*sc2)
945 >             duduz_i(1) = duduz_i(1) + sw*pref*ri3 &
946 >                                         *(uz_j(1) - 3.0d0*ct_j*xhat*sc2)
947 >             duduz_i(2) = duduz_i(2) + sw*pref*ri3 &
948 >                                         *(uz_j(2) - 3.0d0*ct_j*yhat*sc2)
949 >             duduz_i(3) = duduz_i(3) + sw*pref*ri3 &
950 >                                         *(uz_j(3) - 3.0d0*ct_j*zhat*sc2)
951              
952 <             duduz_j(1) = duduz_j(1) + pref*ri3 &
953 <                  *(uz_i(1) - 3.0d0*ct_i*xhat*sc2)
954 <             duduz_j(2) = duduz_j(2) + pref*ri3 &
955 <                  *(uz_i(2) - 3.0d0*ct_i*yhat*sc2)
956 <             duduz_j(3) = duduz_j(3) + pref*ri3 &
957 <                  *(uz_i(3) - 3.0d0*ct_i*zhat*sc2)
952 >             duduz_j(1) = duduz_j(1) + sw*pref*ri3 &
953 >                                         *(uz_i(1) - 3.0d0*ct_i*xhat*sc2)
954 >             duduz_j(2) = duduz_j(2) + sw*pref*ri3 &
955 >                                         *(uz_i(2) - 3.0d0*ct_i*yhat*sc2)
956 >             duduz_j(3) = duduz_j(3) + sw*pref*ri3 &
957 >                                         *(uz_i(3) - 3.0d0*ct_i*zhat*sc2)
958            endif
959         endif
960      endif
# Line 956 | Line 969 | contains
969            cy2 = cy_i * cy_i
970            cz2 = cz_i * cz_i
971  
972 <          pref = sw * pre14 * q_j / 3.0_dp
973 <
961 <          if (corrMethod .eq. 1) then
972 >          if (summationMethod .eq. UNDAMPED_WOLF) then
973 >             pref = pre14 * q_j / 3.0_dp
974               vterm1 = pref * ri3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
975                    qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
976                    qzz_i * (3.0_dp*cz2 - 1.0_dp) )
977               vterm2 = pref * rcuti3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
978                    qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
979                    qzz_i * (3.0_dp*cz2 - 1.0_dp) )
980 <             vpair = vpair + swi * ( vterm1 - vterm2 )
981 <             epot = epot + ( vterm1 - vterm2 )
980 >             vpair = vpair + ( vterm1 - vterm2 )
981 >             epot = epot + sw*( vterm1 - vterm2 )
982              
983 <             dudx = dudx - (5.0_dp*(vterm1*riji*xhat - vterm2*rcuti2*d(1))) + &
984 <                  pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(1)) - &
983 >             dudx = dudx - sw*(5.0_dp*(vterm1*riji*xhat-vterm2*rcuti2*d(1))) +&
984 >                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(1)) - &
985                    qxx_i*2.0_dp*(xhat - rcuti*d(1))) + &
986                    (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(1)) - &
987                    qyy_i*2.0_dp*(xhat - rcuti*d(1))) + &
988                    (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(1)) - &
989                    qzz_i*2.0_dp*(xhat - rcuti*d(1))) )
990 <             dudy = dudy - (5.0_dp*(vterm1*riji*yhat - vterm2*rcuti2*d(2))) + &
991 <                  pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(2)) - &
990 >             dudy = dudy - sw*(5.0_dp*(vterm1*riji*yhat-vterm2*rcuti2*d(2))) +&
991 >                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(2)) - &
992                    qxx_i*2.0_dp*(yhat - rcuti*d(2))) + &
993                    (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(2)) - &
994                    qyy_i*2.0_dp*(yhat - rcuti*d(2))) + &
995                    (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(2)) - &
996                    qzz_i*2.0_dp*(yhat - rcuti*d(2))) )
997 <             dudz = dudz - (5.0_dp*(vterm1*riji*zhat - vterm2*rcuti2*d(3))) + &
998 <                  pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(3)) - &
997 >             dudz = dudz - sw*(5.0_dp*(vterm1*riji*zhat-vterm2*rcuti2*d(3))) +&
998 >                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(3)) - &
999                    qxx_i*2.0_dp*(zhat - rcuti*d(3))) + &
1000                    (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(3)) - &
1001                    qyy_i*2.0_dp*(zhat - rcuti*d(3))) + &
1002                    (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(3)) - &
1003                    qzz_i*2.0_dp*(zhat - rcuti*d(3))) )
1004              
1005 <             dudux_i(1) = dudux_i(1) + pref * (ri3*(qxx_i*6.0_dp*cx_i*xhat) - &
1005 >             dudux_i(1) = dudux_i(1) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*xhat) -&
1006                    rcuti4*(qxx_i*6.0_dp*cx_i*d(1)))
1007 <             dudux_i(2) = dudux_i(2) + pref * (ri3*(qxx_i*6.0_dp*cx_i*yhat) - &
1007 >             dudux_i(2) = dudux_i(2) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*yhat) -&
1008                    rcuti4*(qxx_i*6.0_dp*cx_i*d(2)))
1009 <             dudux_i(3) = dudux_i(3) + pref * (ri3*(qxx_i*6.0_dp*cx_i*zhat) - &
1009 >             dudux_i(3) = dudux_i(3) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*zhat) -&
1010                    rcuti4*(qxx_i*6.0_dp*cx_i*d(3)))
1011              
1012 <             duduy_i(1) = duduy_i(1) + pref * (ri3*(qyy_i*6.0_dp*cy_i*xhat) - &
1012 >             duduy_i(1) = duduy_i(1) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*xhat) -&
1013                    rcuti4*(qyy_i*6.0_dp*cx_i*d(1)))
1014 <             duduy_i(2) = duduy_i(2) + pref * (ri3*(qyy_i*6.0_dp*cy_i*yhat) - &
1014 >             duduy_i(2) = duduy_i(2) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*yhat) -&
1015                    rcuti4*(qyy_i*6.0_dp*cx_i*d(2)))
1016 <             duduy_i(3) = duduy_i(3) + pref * (ri3*(qyy_i*6.0_dp*cy_i*zhat) - &
1016 >             duduy_i(3) = duduy_i(3) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*zhat) -&
1017                    rcuti4*(qyy_i*6.0_dp*cx_i*d(3)))
1018              
1019 <             duduz_i(1) = duduz_i(1) + pref * (ri3*(qzz_i*6.0_dp*cz_i*xhat) - &
1019 >             duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*xhat) -&
1020                    rcuti4*(qzz_i*6.0_dp*cx_i*d(1)))
1021 <             duduz_i(2) = duduz_i(2) + pref * (ri3*(qzz_i*6.0_dp*cz_i*yhat) - &
1021 >             duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*yhat) -&
1022                    rcuti4*(qzz_i*6.0_dp*cx_i*d(2)))
1023 <             duduz_i(3) = duduz_i(3) + pref * (ri3*(qzz_i*6.0_dp*cz_i*zhat) - &
1023 >             duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*zhat) -&
1024                    rcuti4*(qzz_i*6.0_dp*cx_i*d(3)))
1025  
1026            else
1027 +             pref = pre14 * q_j / 3.0_dp
1028               vterm = pref * ri3 * (qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1029                    qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1030                    qzz_i * (3.0_dp*cz2 - 1.0_dp))
1031 <             vpair = vpair + swi * vterm
1032 <             epot = epot + vterm
1031 >             vpair = vpair + vterm
1032 >             epot = epot + sw*vterm
1033              
1034 <             dudx = dudx - 5.0_dp*vterm*riji*xhat + pref * ri4 * ( &
1034 >             dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref*ri4 * ( &
1035                    qxx_i*(6.0_dp*cx_i*ux_i(1) - 2.0_dp*xhat) + &
1036                    qyy_i*(6.0_dp*cy_i*uy_i(1) - 2.0_dp*xhat) + &
1037                    qzz_i*(6.0_dp*cz_i*uz_i(1) - 2.0_dp*xhat) )
1038 <             dudy = dudy - 5.0_dp*vterm*riji*yhat + pref * ri4 * ( &
1038 >             dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref*ri4 * ( &
1039                    qxx_i*(6.0_dp*cx_i*ux_i(2) - 2.0_dp*yhat) + &
1040                    qyy_i*(6.0_dp*cy_i*uy_i(2) - 2.0_dp*yhat) + &
1041                    qzz_i*(6.0_dp*cz_i*uz_i(2) - 2.0_dp*yhat) )
1042 <             dudz = dudz - 5.0_dp*vterm*riji*zhat + pref * ri4 * ( &
1042 >             dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref*ri4 * ( &
1043                    qxx_i*(6.0_dp*cx_i*ux_i(3) - 2.0_dp*zhat) + &
1044                    qyy_i*(6.0_dp*cy_i*uy_i(3) - 2.0_dp*zhat) + &
1045                    qzz_i*(6.0_dp*cz_i*uz_i(3) - 2.0_dp*zhat) )
1046              
1047 <             dudux_i(1) = dudux_i(1) + pref * ri3*(qxx_i*6.0_dp*cx_i*xhat)
1048 <             dudux_i(2) = dudux_i(2) + pref * ri3*(qxx_i*6.0_dp*cx_i*yhat)
1049 <             dudux_i(3) = dudux_i(3) + pref * ri3*(qxx_i*6.0_dp*cx_i*zhat)
1047 >             dudux_i(1) = dudux_i(1) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*xhat)
1048 >             dudux_i(2) = dudux_i(2) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*yhat)
1049 >             dudux_i(3) = dudux_i(3) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*zhat)
1050              
1051 <             duduy_i(1) = duduy_i(1) + pref * ri3*(qyy_i*6.0_dp*cy_i*xhat)
1052 <             duduy_i(2) = duduy_i(2) + pref * ri3*(qyy_i*6.0_dp*cy_i*yhat)
1053 <             duduy_i(3) = duduy_i(3) + pref * ri3*(qyy_i*6.0_dp*cy_i*zhat)
1051 >             duduy_i(1) = duduy_i(1) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*xhat)
1052 >             duduy_i(2) = duduy_i(2) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*yhat)
1053 >             duduy_i(3) = duduy_i(3) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*zhat)
1054              
1055 <             duduz_i(1) = duduz_i(1) + pref * ri3*(qzz_i*6.0_dp*cz_i*xhat)
1056 <             duduz_i(2) = duduz_i(2) + pref * ri3*(qzz_i*6.0_dp*cz_i*yhat)
1057 <             duduz_i(3) = duduz_i(3) + pref * ri3*(qzz_i*6.0_dp*cz_i*zhat)
1055 >             duduz_i(1) = duduz_i(1) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*xhat)
1056 >             duduz_i(2) = duduz_i(2) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*yhat)
1057 >             duduz_i(3) = duduz_i(3) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*zhat)
1058            endif
1059         endif
1060      endif

Diff Legend

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