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 2325 by chrisfen, Sat Sep 24 17:39:36 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  
94    public :: setElectrostaticSummationMethod
# Line 118 | Line 124 | contains
124   contains
125  
126    subroutine setElectrostaticSummationMethod(the_ESM)
121 –
127      integer, intent(in) :: the_ESM    
128  
129      if ((the_ESM .le. 0) .or. (the_ESM .gt. REACTION_FIELD)) then
130         call handleError("setElectrostaticSummationMethod", "Unsupported Summation Method")
131      endif
132 +
133 +    summationMethod = the_ESM
134  
135 +    if (summationMethod .eq. UNDAMPED_WOLF) is_Undamped_Wolf = .true.
136 +    if (summationMethod .eq. DAMPED_WOLF) is_Damped_Wolf = .true.
137    end subroutine setElectrostaticSummationMethod
138  
139    subroutine setElectrostaticCutoffRadius(thisRcut)
# Line 357 | Line 366 | contains
366  
367    subroutine checkSummationMethod()
368  
369 +    if (.not.haveDefaultCutoff) then
370 +       call handleError("checkSummationMethod", "no Default Cutoff set!")
371 +    endif
372 +
373 +    rcuti = 1.0d0 / defaultCutoff
374 +    rcuti2 = rcuti*rcuti
375 +    rcuti3 = rcuti2*rcuti
376 +    rcuti4 = rcuti2*rcuti2
377 +
378      if (summationMethod .eq. DAMPED_WOLF) then
379         if (.not.haveDWAconstants) then
380            
# Line 387 | Line 405 | contains
405  
406  
407    subroutine doElectrostaticPair(atom1, atom2, d, rij, r2, sw, &
408 <       vpair, fpair, pot, eFrame, f, t, do_pot, corrMethod, rcuti)
408 >       vpair, fpair, pot, eFrame, f, t, do_pot)
409  
410      logical, intent(in) :: do_pot
411  
412      integer, intent(in) :: atom1, atom2
413      integer :: localError
396 –    integer, intent(in) :: corrMethod
414  
415 <    real(kind=dp), intent(in) :: rij, r2, sw, rcuti
415 >    real(kind=dp), intent(in) :: rij, r2, sw
416      real(kind=dp), intent(in), dimension(3) :: d
417      real(kind=dp), intent(inout) :: vpair
418      real(kind=dp), intent(inout), dimension(3) :: fpair
419  
420 <    real( kind = dp ) :: pot, swi
420 >    real( kind = dp ) :: pot
421      real( kind = dp ), dimension(9,nLocal) :: eFrame
422      real( kind = dp ), dimension(3,nLocal) :: f
423      real( kind = dp ), dimension(3,nLocal) :: t
# Line 425 | Line 442 | contains
442      real (kind=dp) :: pref, vterm, epot, dudr, vterm1, vterm2
443      real (kind=dp) :: xhat, yhat, zhat
444      real (kind=dp) :: dudx, dudy, dudz
445 <    real (kind=dp) :: scale, sc2, bigR, switcher, dswitcher
429 <    real (kind=dp) :: rcuti2, rcuti3, rcuti4
445 >    real (kind=dp) :: scale, sc2, bigR
446  
447      if (.not.allocated(ElectrostaticMap)) then
448         call handleError("electrostatic", "no ElectrostaticMap was present before first call of do_electrostatic_pair!")
# Line 435 | Line 451 | contains
451  
452      if (.not.summationMethodChecked) then
453         call checkSummationMethod()
454 +      
455      endif
456  
457  
# Line 454 | Line 471 | contains
471      yhat = d(2) * riji
472      zhat = d(3) * riji
473  
457 –    rcuti2 = rcuti*rcuti
458 –    rcuti3 = rcuti2*rcuti
459 –    rcuti4 = rcuti2*rcuti2
460 –
461 –    swi = 1.0d0 / sw
462 –
474      !! logicals
475      i_is_Charge = ElectrostaticMap(me1)%is_Charge
476      i_is_Dipole = ElectrostaticMap(me1)%is_Dipole
# Line 578 | Line 589 | contains
589         cz_j = uz_j(1)*xhat + uz_j(2)*yhat + uz_j(3)*zhat
590      endif
591    
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 –
592      epot = 0.0_dp
593      dudx = 0.0_dp
594      dudy = 0.0_dp
# Line 603 | Line 606 | contains
606  
607         if (j_is_Charge) then
608  
609 <          if (corrMethod .eq. 1) then
607 <             vterm = pre11 * q_i * q_j * (riji - rcuti)
609 >          if (summationMethod .eq. UNDAMPED_WOLF) then
610  
611 +             vterm = pre11 * q_i * q_j * (riji - rcuti)
612               vpair = vpair + vterm
613 <             epot = epot + sw * vterm
613 >             epot = epot + sw*vterm
614              
615               dudr  = - sw * pre11 * q_i * q_j * (riji*riji*riji - rcuti2*rcuti)
616              
# Line 616 | Line 619 | contains
619               dudz = dudz + dudr * d(3)
620  
621            else
619 –             vterm = pre11 * q_i * q_j * riji
622  
623 +             vterm = pre11 * q_i * q_j * riji
624               vpair = vpair + vterm
625 <             epot = epot + sw * vterm
625 >             epot = epot + sw*vterm
626              
627               dudr  = - sw * vterm * riji
628              
# Line 633 | Line 636 | contains
636  
637         if (j_is_Dipole) then
638  
639 <          pref = sw * pre12 * q_i * mu_j
639 >          pref = pre12 * q_i * mu_j
640  
641 <          if (corrMethod .eq. 1) then
641 >          if (summationMethod .eq. UNDAMPED_WOLF) then
642               ri2 = riji * riji
643               ri3 = ri2 * riji
644  
645 +             pref = pre12 * q_i * mu_j
646               vterm = - pref * ct_j * (ri2 - rcuti2)
647 <             vpair = vpair + swi*vterm
648 <             epot = epot + vterm
647 >             vpair = vpair + vterm
648 >             epot = epot + sw*vterm
649              
650               !! this has a + sign in the () because the rij vector is
651               !! r_j - r_i and the charge-dipole potential takes the origin
652               !! as the point dipole, which is atom j in this case.
653              
654 <             dudx = dudx - pref * ( ri3*( uz_j(1) - 3.0d0*ct_j*xhat) &
654 >             dudx = dudx - sw*pref * ( ri3*( uz_j(1) - 3.0d0*ct_j*xhat) &
655                    - rcuti3*( uz_j(1) - 3.0d0*ct_j*d(1)*rcuti ) )
656 <             dudy = dudy - pref * ( ri3*( uz_j(2) - 3.0d0*ct_j*yhat) &
656 >             dudy = dudy - sw*pref * ( ri3*( uz_j(2) - 3.0d0*ct_j*yhat) &
657                    - rcuti3*( uz_j(2) - 3.0d0*ct_j*d(2)*rcuti ) )
658 <             dudz = dudz - pref * ( ri3*( uz_j(3) - 3.0d0*ct_j*zhat) &
658 >             dudz = dudz - sw*pref * ( ri3*( uz_j(3) - 3.0d0*ct_j*zhat) &
659                    - rcuti3*( uz_j(3) - 3.0d0*ct_j*d(3)*rcuti ) )
660              
661 <             duduz_j(1) = duduz_j(1) - pref*( ri2*xhat - d(1)*rcuti3 )
662 <             duduz_j(2) = duduz_j(2) - pref*( ri2*yhat - d(2)*rcuti3 )
663 <             duduz_j(3) = duduz_j(3) - pref*( ri2*zhat - d(3)*rcuti3 )
661 >             duduz_j(1) = duduz_j(1) - sw*pref*( ri2*xhat - d(1)*rcuti3 )
662 >             duduz_j(2) = duduz_j(2) - sw*pref*( ri2*yhat - d(2)*rcuti3 )
663 >             duduz_j(3) = duduz_j(3) - sw*pref*( ri2*zhat - d(3)*rcuti3 )
664  
665            else
666               if (j_is_SplitDipole) then
# Line 671 | Line 675 | contains
675               ri2 = ri * ri
676               ri3 = ri2 * ri
677               sc2 = scale * scale
678 <            
678 >
679 >             pref = pre12 * q_i * mu_j
680               vterm = - pref * ct_j * ri2 * scale
681 <             vpair = vpair + swi * vterm
682 <             epot = epot + vterm
681 >             vpair = vpair + vterm
682 >             epot = epot + sw*vterm
683              
684               !! this has a + sign in the () because the rij vector is
685               !! r_j - r_i and the charge-dipole potential takes the origin
686               !! as the point dipole, which is atom j in this case.
687              
688 <             dudx = dudx - pref * ri3 * ( uz_j(1) - 3.0d0*ct_j*xhat*sc2)
689 <             dudy = dudy - pref * ri3 * ( uz_j(2) - 3.0d0*ct_j*yhat*sc2)
690 <             dudz = dudz - pref * ri3 * ( uz_j(3) - 3.0d0*ct_j*zhat*sc2)
688 >             dudx = dudx - sw*pref * ri3 * ( uz_j(1) - 3.0d0*ct_j*xhat*sc2)
689 >             dudy = dudy - sw*pref * ri3 * ( uz_j(2) - 3.0d0*ct_j*yhat*sc2)
690 >             dudz = dudz - sw*pref * ri3 * ( uz_j(3) - 3.0d0*ct_j*zhat*sc2)
691              
692 <             duduz_j(1) = duduz_j(1) - pref * ri2 * xhat * scale
693 <             duduz_j(2) = duduz_j(2) - pref * ri2 * yhat * scale
694 <             duduz_j(3) = duduz_j(3) - pref * ri2 * zhat * scale
692 >             duduz_j(1) = duduz_j(1) - sw*pref * ri2 * xhat * scale
693 >             duduz_j(2) = duduz_j(2) - sw*pref * ri2 * yhat * scale
694 >             duduz_j(3) = duduz_j(3) - sw*pref * ri2 * zhat * scale
695  
696            endif
697         endif
# Line 699 | Line 704 | contains
704            cy2 = cy_j * cy_j
705            cz2 = cz_j * cz_j
706  
707 <
708 <          pref =  sw * pre14 * q_i / 3.0_dp
704 <
705 <          if (corrMethod .eq. 1) then
707 >          if (summationMethod .eq. UNDAMPED_WOLF) then
708 >             pref =  pre14 * q_i / 3.0_dp
709               vterm1 = pref * ri3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
710                    qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
711                    qzz_j * (3.0_dp*cz2 - 1.0_dp) )
712               vterm2 = pref * rcuti3*( 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 <             vpair = vpair + swi*( vterm1 - vterm2 )
716 <             epot = epot + ( vterm1 - vterm2 )
715 >             vpair = vpair + ( vterm1 - vterm2 )
716 >             epot = epot + sw*( vterm1 - vterm2 )
717              
718               dudx = dudx - (5.0_dp * &
719 <                  (vterm1*riji*xhat - vterm2*rcuti2*d(1))) + pref * ( &
719 >                  (vterm1*riji*xhat - vterm2*rcuti2*d(1))) + sw*pref * ( &
720                    (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(1)) - &
721                    qxx_j*2.0_dp*(xhat - rcuti*d(1))) + &
722                    (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(1)) - &
# Line 721 | Line 724 | contains
724                    (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(1)) - &
725                    qzz_j*2.0_dp*(xhat - rcuti*d(1))) )
726               dudy = dudy - (5.0_dp * &
727 <                  (vterm1*riji*yhat - vterm2*rcuti2*d(2))) + pref * ( &
727 >                  (vterm1*riji*yhat - vterm2*rcuti2*d(2))) + sw*pref * ( &
728                    (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(2)) - &
729                    qxx_j*2.0_dp*(yhat - rcuti*d(2))) + &
730                    (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(2)) - &
# Line 729 | Line 732 | contains
732                    (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(2)) - &
733                    qzz_j*2.0_dp*(yhat - rcuti*d(2))) )
734               dudz = dudz - (5.0_dp * &
735 <                  (vterm1*riji*zhat - vterm2*rcuti2*d(3))) + pref * ( &
735 >                  (vterm1*riji*zhat - vterm2*rcuti2*d(3))) + sw*pref * ( &
736                    (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(3)) - &
737                    qxx_j*2.0_dp*(zhat - rcuti*d(3))) + &
738                    (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(3)) - &
# Line 737 | Line 740 | contains
740                    (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(3)) - &
741                    qzz_j*2.0_dp*(zhat - rcuti*d(3))) )
742              
743 <             dudux_j(1) = dudux_j(1) + pref * (ri3*(qxx_j*6.0_dp*cx_j*xhat) - &
743 >             dudux_j(1) = dudux_j(1) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*xhat) -&
744                    rcuti4*(qxx_j*6.0_dp*cx_j*d(1)))
745 <             dudux_j(2) = dudux_j(2) + pref * (ri3*(qxx_j*6.0_dp*cx_j*yhat) - &
745 >             dudux_j(2) = dudux_j(2) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*yhat) -&
746                    rcuti4*(qxx_j*6.0_dp*cx_j*d(2)))
747 <             dudux_j(3) = dudux_j(3) + pref * (ri3*(qxx_j*6.0_dp*cx_j*zhat) - &
747 >             dudux_j(3) = dudux_j(3) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*zhat) -&
748                    rcuti4*(qxx_j*6.0_dp*cx_j*d(3)))
749              
750 <             duduy_j(1) = duduy_j(1) + pref * (ri3*(qyy_j*6.0_dp*cy_j*xhat) - &
750 >             duduy_j(1) = duduy_j(1) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*xhat) -&
751                    rcuti4*(qyy_j*6.0_dp*cx_j*d(1)))
752 <             duduy_j(2) = duduy_j(2) + pref * (ri3*(qyy_j*6.0_dp*cy_j*yhat) - &
752 >             duduy_j(2) = duduy_j(2) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*yhat) -&
753                    rcuti4*(qyy_j*6.0_dp*cx_j*d(2)))
754 <             duduy_j(3) = duduy_j(3) + pref * (ri3*(qyy_j*6.0_dp*cy_j*zhat) - &
754 >             duduy_j(3) = duduy_j(3) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*zhat) -&
755                    rcuti4*(qyy_j*6.0_dp*cx_j*d(3)))
756              
757 <             duduz_j(1) = duduz_j(1) + pref * (ri3*(qzz_j*6.0_dp*cz_j*xhat) - &
757 >             duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*xhat) -&
758                    rcuti4*(qzz_j*6.0_dp*cx_j*d(1)))
759 <             duduz_j(2) = duduz_j(2) + pref * (ri3*(qzz_j*6.0_dp*cz_j*yhat) - &
759 >             duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*yhat) -&
760                    rcuti4*(qzz_j*6.0_dp*cx_j*d(2)))
761 <             duduz_j(3) = duduz_j(3) + pref * (ri3*(qzz_j*6.0_dp*cz_j*zhat) - &
761 >             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*zhat) -&
762                    rcuti4*(qzz_j*6.0_dp*cx_j*d(3)))
763          
764            else
765 +             pref =  pre14 * q_i / 3.0_dp
766               vterm = pref * ri3 * (qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
767                    qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
768                    qzz_j * (3.0_dp*cz2 - 1.0_dp))
769 <             vpair = vpair + swi * vterm
770 <             epot = epot + vterm
769 >             vpair = vpair + vterm
770 >             epot = epot + sw*vterm
771              
772 <             dudx = dudx - 5.0_dp*vterm*riji*xhat + pref * ri4 * ( &
772 >             dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref * ri4 * ( &
773                    qxx_j*(6.0_dp*cx_j*ux_j(1) - 2.0_dp*xhat) + &
774                    qyy_j*(6.0_dp*cy_j*uy_j(1) - 2.0_dp*xhat) + &
775                    qzz_j*(6.0_dp*cz_j*uz_j(1) - 2.0_dp*xhat) )
776 <             dudy = dudy - 5.0_dp*vterm*riji*yhat + pref * ri4 * ( &
776 >             dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref * ri4 * ( &
777                    qxx_j*(6.0_dp*cx_j*ux_j(2) - 2.0_dp*yhat) + &
778                    qyy_j*(6.0_dp*cy_j*uy_j(2) - 2.0_dp*yhat) + &
779                    qzz_j*(6.0_dp*cz_j*uz_j(2) - 2.0_dp*yhat) )
780 <             dudz = dudz - 5.0_dp*vterm*riji*zhat + pref * ri4 * ( &
780 >             dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref * ri4 * ( &
781                    qxx_j*(6.0_dp*cx_j*ux_j(3) - 2.0_dp*zhat) + &
782                    qyy_j*(6.0_dp*cy_j*uy_j(3) - 2.0_dp*zhat) + &
783                    qzz_j*(6.0_dp*cz_j*uz_j(3) - 2.0_dp*zhat) )
784              
785 <             dudux_j(1) = dudux_j(1) + pref * ri3*(qxx_j*6.0_dp*cx_j*xhat)
786 <             dudux_j(2) = dudux_j(2) + pref * ri3*(qxx_j*6.0_dp*cx_j*yhat)
787 <             dudux_j(3) = dudux_j(3) + pref * ri3*(qxx_j*6.0_dp*cx_j*zhat)
785 >             dudux_j(1) = dudux_j(1) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*xhat)
786 >             dudux_j(2) = dudux_j(2) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*yhat)
787 >             dudux_j(3) = dudux_j(3) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*zhat)
788              
789 <             duduy_j(1) = duduy_j(1) + pref * ri3*(qyy_j*6.0_dp*cy_j*xhat)
790 <             duduy_j(2) = duduy_j(2) + pref * ri3*(qyy_j*6.0_dp*cy_j*yhat)
791 <             duduy_j(3) = duduy_j(3) + pref * ri3*(qyy_j*6.0_dp*cy_j*zhat)
789 >             duduy_j(1) = duduy_j(1) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*xhat)
790 >             duduy_j(2) = duduy_j(2) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*yhat)
791 >             duduy_j(3) = duduy_j(3) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*zhat)
792              
793 <             duduz_j(1) = duduz_j(1) + pref * ri3*(qzz_j*6.0_dp*cz_j*xhat)
794 <             duduz_j(2) = duduz_j(2) + pref * ri3*(qzz_j*6.0_dp*cz_j*yhat)
795 <             duduz_j(3) = duduz_j(3) + pref * ri3*(qzz_j*6.0_dp*cz_j*zhat)
793 >             duduz_j(1) = duduz_j(1) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*xhat)
794 >             duduz_j(2) = duduz_j(2) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*yhat)
795 >             duduz_j(3) = duduz_j(3) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*zhat)
796            
797            endif
798         endif
# Line 797 | Line 801 | contains
801      if (i_is_Dipole) then
802  
803         if (j_is_Charge) then
804 <
805 <          pref = sw * pre12 * q_j * mu_i
806 <
807 <          if (corrMethod .eq. 1) then
804 >          
805 >          pref = pre12 * q_j * mu_i
806 >          
807 >          if (summationMethod .eq. UNDAMPED_WOLF) then
808               ri2 = riji * riji
809               ri3 = ri2 * riji
810  
811 +             pref = pre12 * q_j * mu_i
812               vterm = pref * ct_i * (ri2 - rcuti2)
813 <             vpair = vpair + swi * vterm
814 <             epot = epot + vterm
813 >             vpair = vpair + vterm
814 >             epot = epot + sw*vterm
815              
816               !! this has a + sign in the () because the rij vector is
817               !! r_j - r_i and the charge-dipole potential takes the origin
818               !! as the point dipole, which is atom j in this case.
819              
820 <             dudx = dudx + pref * ( ri3*( uz_i(1) - 3.0d0*ct_i*xhat) &
820 >             dudx = dudx + sw*pref * ( ri3*( uz_i(1) - 3.0d0*ct_i*xhat) &
821                    - rcuti3*( uz_i(1) - 3.0d0*ct_i*d(1)*rcuti ) )
822 <             dudy = dudy + pref * ( ri3*( uz_i(2) - 3.0d0*ct_i*yhat) &
822 >             dudy = dudy + sw*pref * ( ri3*( uz_i(2) - 3.0d0*ct_i*yhat) &
823                    - rcuti3*( uz_i(2) - 3.0d0*ct_i*d(2)*rcuti ) )
824 <             dudz = dudz + pref * ( ri3*( uz_i(3) - 3.0d0*ct_i*zhat) &
824 >             dudz = dudz + sw*pref * ( ri3*( uz_i(3) - 3.0d0*ct_i*zhat) &
825                    - rcuti3*( uz_i(3) - 3.0d0*ct_i*d(3)*rcuti ) )
826              
827 <             duduz_i(1) = duduz_i(1) - pref*( ri2*xhat - d(1)*rcuti3 )
828 <             duduz_i(2) = duduz_i(2) - pref*( ri2*yhat - d(2)*rcuti3 )
829 <             duduz_i(3) = duduz_i(3) - pref*( ri2*zhat - d(3)*rcuti3 )
827 >             duduz_i(1) = duduz_i(1) - sw*pref*( ri2*xhat - d(1)*rcuti3 )
828 >             duduz_i(2) = duduz_i(2) - sw*pref*( ri2*yhat - d(2)*rcuti3 )
829 >             duduz_i(3) = duduz_i(3) - sw*pref*( ri2*zhat - d(3)*rcuti3 )
830  
831            else
832               if (i_is_SplitDipole) then
# Line 836 | Line 841 | contains
841               ri2 = ri * ri
842               ri3 = ri2 * ri
843               sc2 = scale * scale
844 <            
844 >
845 >             pref = pre12 * q_j * mu_i
846               vterm = pref * ct_i * ri2 * scale
847 <             vpair = vpair + swi * vterm
848 <             epot = epot + vterm
847 >             vpair = vpair + vterm
848 >             epot = epot + sw*vterm
849              
850 <             dudx = dudx + pref * ri3 * ( uz_i(1) - 3.0d0 * ct_i * xhat*sc2)
851 <             dudy = dudy + pref * ri3 * ( uz_i(2) - 3.0d0 * ct_i * yhat*sc2)
852 <             dudz = dudz + pref * ri3 * ( uz_i(3) - 3.0d0 * ct_i * zhat*sc2)
850 >             dudx = dudx + sw*pref * ri3 * ( uz_i(1) - 3.0d0 * ct_i * xhat*sc2)
851 >             dudy = dudy + sw*pref * ri3 * ( uz_i(2) - 3.0d0 * ct_i * yhat*sc2)
852 >             dudz = dudz + sw*pref * ri3 * ( uz_i(3) - 3.0d0 * ct_i * zhat*sc2)
853              
854 <             duduz_i(1) = duduz_i(1) + pref * ri2 * xhat * scale
855 <             duduz_i(2) = duduz_i(2) + pref * ri2 * yhat * scale
856 <             duduz_i(3) = duduz_i(3) + pref * ri2 * zhat * scale
854 >             duduz_i(1) = duduz_i(1) + sw*pref * ri2 * xhat * scale
855 >             duduz_i(2) = duduz_i(2) + sw*pref * ri2 * yhat * scale
856 >             duduz_i(3) = duduz_i(3) + sw*pref * ri2 * zhat * scale
857            endif
858         endif
859 <
859 >      
860         if (j_is_Dipole) then
861  
862 <          pref = sw * pre22 * mu_i * mu_j
857 <
858 <          if (corrMethod .eq. 1) then
862 >          if (summationMethod .eq. UNDAMPED_WOLF) then
863               ri2 = riji * riji
864               ri3 = ri2 * riji
865               ri4 = ri2 * ri2
866  
867 +             pref = pre22 * mu_i * mu_j
868               vterm = pref * (ri3 - rcuti3) * (ct_ij - 3.0d0 * ct_i * ct_j)
869 <             vpair = vpair + swi * vterm
870 <             epot = epot + vterm
869 >             vpair = vpair + vterm
870 >             epot = epot + sw*vterm
871              
872               a1 = 5.0d0 * ct_i * ct_j - ct_ij
873              
874 <             dudx = dudx + pref*3.0d0*ri4 &
875 <                  *(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1)) - &
876 <                  pref*3.0d0*rcuti4*(a1*rcuti*d(1)-ct_i*uz_j(1)-ct_j*uz_i(1))
877 <             dudy = dudy + pref*3.0d0*ri4 &
878 <                  *(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2)) - &
879 <                  pref*3.0d0*rcuti4*(a1*rcuti*d(2)-ct_i*uz_j(2)-ct_j*uz_i(2))
880 <             dudz = dudz + pref*3.0d0*ri4 &
881 <                  *(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3)) - &
882 <                  pref*3.0d0*rcuti4*(a1*rcuti*d(3)-ct_i*uz_j(3)-ct_j*uz_i(3))
874 >             dudx = dudx + sw*pref*3.0d0*ri4 &
875 >                             * (a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1)) &
876 >                         - sw*pref*3.0d0*rcuti4 &
877 >                             * (a1*rcuti*d(1)-ct_i*uz_j(1)-ct_j*uz_i(1))
878 >             dudy = dudy + sw*pref*3.0d0*ri4 &
879 >                             * (a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2)) &
880 >                         - sw*pref*3.0d0*rcuti4 &
881 >                             * (a1*rcuti*d(2)-ct_i*uz_j(2)-ct_j*uz_i(2))
882 >             dudz = dudz + sw*pref*3.0d0*ri4 &
883 >                             * (a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3)) &
884 >                         - sw*pref*3.0d0*rcuti4 &
885 >                             * (a1*rcuti*d(3)-ct_i*uz_j(3)-ct_j*uz_i(3))
886              
887 <             duduz_i(1) = duduz_i(1) + pref*(ri3*(uz_j(1) - 3.0d0*ct_j*xhat) &
887 >             duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(uz_j(1)-3.0d0*ct_j*xhat) &
888                    - rcuti3*(uz_j(1) - 3.0d0*ct_j*d(1)*rcuti))
889 <             duduz_i(2) = duduz_i(2) + pref*(ri3*(uz_j(2) - 3.0d0*ct_j*yhat) &
889 >             duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(uz_j(2)-3.0d0*ct_j*yhat) &
890                    - rcuti3*(uz_j(2) - 3.0d0*ct_j*d(2)*rcuti))
891 <             duduz_i(3) = duduz_i(3) + pref*(ri3*(uz_j(3) - 3.0d0*ct_j*zhat) &
891 >             duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(uz_j(3)-3.0d0*ct_j*zhat) &
892                    - rcuti3*(uz_j(3) - 3.0d0*ct_j*d(3)*rcuti))
893 <             duduz_j(1) = duduz_j(1) + pref*(ri3*(uz_i(1) - 3.0d0*ct_i*xhat) &
893 >             duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(uz_i(1)-3.0d0*ct_i*xhat) &
894                    - rcuti3*(uz_i(1) - 3.0d0*ct_i*d(1)*rcuti))
895 <             duduz_j(2) = duduz_j(2) + pref*(ri3*(uz_i(2) - 3.0d0*ct_i*yhat) &
895 >             duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(uz_i(2)-3.0d0*ct_i*yhat) &
896                    - rcuti3*(uz_i(2) - 3.0d0*ct_i*d(2)*rcuti))
897 <             duduz_j(3) = duduz_j(3) + pref*(ri3*(uz_i(3) - 3.0d0*ct_i*zhat) &
897 >             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(uz_i(3)-3.0d0*ct_i*zhat) &
898                    - rcuti3*(uz_i(3) - 3.0d0*ct_i*d(3)*rcuti))
899 +
900            else
892 –            
901               if (i_is_SplitDipole) then
902                  if (j_is_SplitDipole) then
903                     BigR = sqrt(r2 + 0.25_dp * d_i * d_i + 0.25_dp * d_j * d_j)
# Line 916 | Line 924 | contains
924               ri4 = ri2 * ri2
925               sc2 = scale * scale
926              
927 +             pref = pre22 * mu_i * mu_j
928               vterm = pref * ri3 * (ct_ij - 3.0d0 * ct_i * ct_j * sc2)
929 <             vpair = vpair + swi * vterm
930 <             epot = epot + vterm
929 >             vpair = vpair + vterm
930 >             epot = epot + sw*vterm
931              
932               a1 = 5.0d0 * ct_i * ct_j * sc2 - ct_ij
933              
934 <             dudx = dudx + pref*3.0d0*ri4*scale &
935 <                  *(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
936 <             dudy = dudy + pref*3.0d0*ri4*scale &
937 <                  *(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
938 <             dudz = dudz + pref*3.0d0*ri4*scale &
939 <                  *(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
934 >             dudx = dudx + sw*pref*3.0d0*ri4*scale &
935 >                             *(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
936 >             dudy = dudy + sw*pref*3.0d0*ri4*scale &
937 >                             *(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
938 >             dudz = dudz + sw*pref*3.0d0*ri4*scale &
939 >                             *(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
940              
941 <             duduz_i(1) = duduz_i(1) + pref*ri3 &
942 <                  *(uz_j(1) - 3.0d0*ct_j*xhat*sc2)
943 <             duduz_i(2) = duduz_i(2) + pref*ri3 &
944 <                  *(uz_j(2) - 3.0d0*ct_j*yhat*sc2)
945 <             duduz_i(3) = duduz_i(3) + pref*ri3 &
946 <                  *(uz_j(3) - 3.0d0*ct_j*zhat*sc2)
941 >             duduz_i(1) = duduz_i(1) + sw*pref*ri3 &
942 >                                         *(uz_j(1) - 3.0d0*ct_j*xhat*sc2)
943 >             duduz_i(2) = duduz_i(2) + sw*pref*ri3 &
944 >                                         *(uz_j(2) - 3.0d0*ct_j*yhat*sc2)
945 >             duduz_i(3) = duduz_i(3) + sw*pref*ri3 &
946 >                                         *(uz_j(3) - 3.0d0*ct_j*zhat*sc2)
947              
948 <             duduz_j(1) = duduz_j(1) + pref*ri3 &
949 <                  *(uz_i(1) - 3.0d0*ct_i*xhat*sc2)
950 <             duduz_j(2) = duduz_j(2) + pref*ri3 &
951 <                  *(uz_i(2) - 3.0d0*ct_i*yhat*sc2)
952 <             duduz_j(3) = duduz_j(3) + pref*ri3 &
953 <                  *(uz_i(3) - 3.0d0*ct_i*zhat*sc2)
948 >             duduz_j(1) = duduz_j(1) + sw*pref*ri3 &
949 >                                         *(uz_i(1) - 3.0d0*ct_i*xhat*sc2)
950 >             duduz_j(2) = duduz_j(2) + sw*pref*ri3 &
951 >                                         *(uz_i(2) - 3.0d0*ct_i*yhat*sc2)
952 >             duduz_j(3) = duduz_j(3) + sw*pref*ri3 &
953 >                                         *(uz_i(3) - 3.0d0*ct_i*zhat*sc2)
954            endif
955         endif
956      endif
# Line 956 | Line 965 | contains
965            cy2 = cy_i * cy_i
966            cz2 = cz_i * cz_i
967  
968 <          pref = sw * pre14 * q_j / 3.0_dp
969 <
961 <          if (corrMethod .eq. 1) then
968 >          if (summationMethod .eq. UNDAMPED_WOLF) then
969 >             pref = pre14 * q_j / 3.0_dp
970               vterm1 = pref * ri3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
971                    qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
972                    qzz_i * (3.0_dp*cz2 - 1.0_dp) )
973               vterm2 = pref * rcuti3*( 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 <             vpair = vpair + swi * ( vterm1 - vterm2 )
977 <             epot = epot + ( vterm1 - vterm2 )
976 >             vpair = vpair + ( vterm1 - vterm2 )
977 >             epot = epot + sw*( vterm1 - vterm2 )
978              
979 <             dudx = dudx - (5.0_dp*(vterm1*riji*xhat - vterm2*rcuti2*d(1))) + &
980 <                  pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(1)) - &
979 >             dudx = dudx - sw*(5.0_dp*(vterm1*riji*xhat-vterm2*rcuti2*d(1))) +&
980 >                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(1)) - &
981                    qxx_i*2.0_dp*(xhat - rcuti*d(1))) + &
982                    (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(1)) - &
983                    qyy_i*2.0_dp*(xhat - rcuti*d(1))) + &
984                    (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(1)) - &
985                    qzz_i*2.0_dp*(xhat - rcuti*d(1))) )
986 <             dudy = dudy - (5.0_dp*(vterm1*riji*yhat - vterm2*rcuti2*d(2))) + &
987 <                  pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(2)) - &
986 >             dudy = dudy - sw*(5.0_dp*(vterm1*riji*yhat-vterm2*rcuti2*d(2))) +&
987 >                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(2)) - &
988                    qxx_i*2.0_dp*(yhat - rcuti*d(2))) + &
989                    (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(2)) - &
990                    qyy_i*2.0_dp*(yhat - rcuti*d(2))) + &
991                    (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(2)) - &
992                    qzz_i*2.0_dp*(yhat - rcuti*d(2))) )
993 <             dudz = dudz - (5.0_dp*(vterm1*riji*zhat - vterm2*rcuti2*d(3))) + &
994 <                  pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(3)) - &
993 >             dudz = dudz - sw*(5.0_dp*(vterm1*riji*zhat-vterm2*rcuti2*d(3))) +&
994 >                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(3)) - &
995                    qxx_i*2.0_dp*(zhat - rcuti*d(3))) + &
996                    (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(3)) - &
997                    qyy_i*2.0_dp*(zhat - rcuti*d(3))) + &
998                    (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(3)) - &
999                    qzz_i*2.0_dp*(zhat - rcuti*d(3))) )
1000              
1001 <             dudux_i(1) = dudux_i(1) + pref * (ri3*(qxx_i*6.0_dp*cx_i*xhat) - &
1001 >             dudux_i(1) = dudux_i(1) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*xhat) -&
1002                    rcuti4*(qxx_i*6.0_dp*cx_i*d(1)))
1003 <             dudux_i(2) = dudux_i(2) + pref * (ri3*(qxx_i*6.0_dp*cx_i*yhat) - &
1003 >             dudux_i(2) = dudux_i(2) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*yhat) -&
1004                    rcuti4*(qxx_i*6.0_dp*cx_i*d(2)))
1005 <             dudux_i(3) = dudux_i(3) + pref * (ri3*(qxx_i*6.0_dp*cx_i*zhat) - &
1005 >             dudux_i(3) = dudux_i(3) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*zhat) -&
1006                    rcuti4*(qxx_i*6.0_dp*cx_i*d(3)))
1007              
1008 <             duduy_i(1) = duduy_i(1) + pref * (ri3*(qyy_i*6.0_dp*cy_i*xhat) - &
1008 >             duduy_i(1) = duduy_i(1) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*xhat) -&
1009                    rcuti4*(qyy_i*6.0_dp*cx_i*d(1)))
1010 <             duduy_i(2) = duduy_i(2) + pref * (ri3*(qyy_i*6.0_dp*cy_i*yhat) - &
1010 >             duduy_i(2) = duduy_i(2) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*yhat) -&
1011                    rcuti4*(qyy_i*6.0_dp*cx_i*d(2)))
1012 <             duduy_i(3) = duduy_i(3) + pref * (ri3*(qyy_i*6.0_dp*cy_i*zhat) - &
1012 >             duduy_i(3) = duduy_i(3) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*zhat) -&
1013                    rcuti4*(qyy_i*6.0_dp*cx_i*d(3)))
1014              
1015 <             duduz_i(1) = duduz_i(1) + pref * (ri3*(qzz_i*6.0_dp*cz_i*xhat) - &
1015 >             duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*xhat) -&
1016                    rcuti4*(qzz_i*6.0_dp*cx_i*d(1)))
1017 <             duduz_i(2) = duduz_i(2) + pref * (ri3*(qzz_i*6.0_dp*cz_i*yhat) - &
1017 >             duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*yhat) -&
1018                    rcuti4*(qzz_i*6.0_dp*cx_i*d(2)))
1019 <             duduz_i(3) = duduz_i(3) + pref * (ri3*(qzz_i*6.0_dp*cz_i*zhat) - &
1019 >             duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*zhat) -&
1020                    rcuti4*(qzz_i*6.0_dp*cx_i*d(3)))
1021  
1022            else
1023 +             pref = pre14 * q_j / 3.0_dp
1024               vterm = pref * ri3 * (qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1025                    qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1026                    qzz_i * (3.0_dp*cz2 - 1.0_dp))
1027 <             vpair = vpair + swi * vterm
1028 <             epot = epot + vterm
1027 >             vpair = vpair + vterm
1028 >             epot = epot + sw*vterm
1029              
1030 <             dudx = dudx - 5.0_dp*vterm*riji*xhat + pref * ri4 * ( &
1030 >             dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref*ri4 * ( &
1031                    qxx_i*(6.0_dp*cx_i*ux_i(1) - 2.0_dp*xhat) + &
1032                    qyy_i*(6.0_dp*cy_i*uy_i(1) - 2.0_dp*xhat) + &
1033                    qzz_i*(6.0_dp*cz_i*uz_i(1) - 2.0_dp*xhat) )
1034 <             dudy = dudy - 5.0_dp*vterm*riji*yhat + pref * ri4 * ( &
1034 >             dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref*ri4 * ( &
1035                    qxx_i*(6.0_dp*cx_i*ux_i(2) - 2.0_dp*yhat) + &
1036                    qyy_i*(6.0_dp*cy_i*uy_i(2) - 2.0_dp*yhat) + &
1037                    qzz_i*(6.0_dp*cz_i*uz_i(2) - 2.0_dp*yhat) )
1038 <             dudz = dudz - 5.0_dp*vterm*riji*zhat + pref * ri4 * ( &
1038 >             dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref*ri4 * ( &
1039                    qxx_i*(6.0_dp*cx_i*ux_i(3) - 2.0_dp*zhat) + &
1040                    qyy_i*(6.0_dp*cy_i*uy_i(3) - 2.0_dp*zhat) + &
1041                    qzz_i*(6.0_dp*cz_i*uz_i(3) - 2.0_dp*zhat) )
1042              
1043 <             dudux_i(1) = dudux_i(1) + pref * ri3*(qxx_i*6.0_dp*cx_i*xhat)
1044 <             dudux_i(2) = dudux_i(2) + pref * ri3*(qxx_i*6.0_dp*cx_i*yhat)
1045 <             dudux_i(3) = dudux_i(3) + pref * ri3*(qxx_i*6.0_dp*cx_i*zhat)
1043 >             dudux_i(1) = dudux_i(1) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*xhat)
1044 >             dudux_i(2) = dudux_i(2) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*yhat)
1045 >             dudux_i(3) = dudux_i(3) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*zhat)
1046              
1047 <             duduy_i(1) = duduy_i(1) + pref * ri3*(qyy_i*6.0_dp*cy_i*xhat)
1048 <             duduy_i(2) = duduy_i(2) + pref * ri3*(qyy_i*6.0_dp*cy_i*yhat)
1049 <             duduy_i(3) = duduy_i(3) + pref * ri3*(qyy_i*6.0_dp*cy_i*zhat)
1047 >             duduy_i(1) = duduy_i(1) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*xhat)
1048 >             duduy_i(2) = duduy_i(2) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*yhat)
1049 >             duduy_i(3) = duduy_i(3) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*zhat)
1050              
1051 <             duduz_i(1) = duduz_i(1) + pref * ri3*(qzz_i*6.0_dp*cz_i*xhat)
1052 <             duduz_i(2) = duduz_i(2) + pref * ri3*(qzz_i*6.0_dp*cz_i*yhat)
1053 <             duduz_i(3) = duduz_i(3) + pref * ri3*(qzz_i*6.0_dp*cz_i*zhat)
1051 >             duduz_i(1) = duduz_i(1) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*xhat)
1052 >             duduz_i(2) = duduz_i(2) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*yhat)
1053 >             duduz_i(3) = duduz_i(3) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*zhat)
1054            endif
1055         endif
1056      endif

Diff Legend

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