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

Comparing trunk/OOPSE-3.0/src/UseTheForce/DarkSide/electrostatic.F90 (file contents):
Revision 2279 by chrisfen, Tue Aug 30 18:23:50 2005 UTC vs.
Revision 2296 by chrisfen, Thu Sep 15 00:13:56 2005 UTC

# Line 307 | Line 307 | contains
307    end function getDipoleMoment
308  
309    subroutine doElectrostaticPair(atom1, atom2, d, rij, r2, sw, &
310 <       vpair, fpair, pot, eFrame, f, t, do_pot, corrMethod)
310 >       vpair, fpair, pot, eFrame, f, t, do_pot, corrMethod, rcuti)
311  
312      logical, intent(in) :: do_pot
313  
# Line 315 | Line 315 | contains
315      integer :: localError
316      integer, intent(in) :: corrMethod
317  
318 <    real(kind=dp), intent(in) :: rij, r2, sw
318 >    real(kind=dp), intent(in) :: rij, r2, sw, rcuti
319      real(kind=dp), intent(in), dimension(3) :: d
320      real(kind=dp), intent(inout) :: vpair
321      real(kind=dp), intent(inout), dimension(3) :: fpair
322  
323 <    real( kind = dp ) :: pot
323 >    real( kind = dp ) :: pot, swi
324      real( kind = dp ), dimension(9,nLocal) :: eFrame
325      real( kind = dp ), dimension(3,nLocal) :: f
326      real( kind = dp ), dimension(3,nLocal) :: t
# Line 333 | Line 333 | contains
333      logical :: i_is_Charge, i_is_Dipole, i_is_SplitDipole, i_is_Quadrupole
334      logical :: j_is_Charge, j_is_Dipole, j_is_SplitDipole, j_is_Quadrupole
335      logical :: i_is_Tap, j_is_Tap
336 –    logical :: use_damped_wolf, use_undamped_wolf
336      integer :: me1, me2, id1, id2
337      real (kind=dp) :: q_i, q_j, mu_i, mu_j, d_i, d_j
338      real (kind=dp) :: qxx_i, qyy_i, qzz_i
# Line 343 | Line 342 | contains
342      real (kind=dp) :: cx2, cy2, cz2
343      real (kind=dp) :: ct_i, ct_j, ct_ij, a1
344      real (kind=dp) :: riji, ri, ri2, ri3, ri4
345 <    real (kind=dp) :: pref, vterm, epot, dudr    
345 >    real (kind=dp) :: pref, vterm, epot, dudr, vterm1, vterm2
346      real (kind=dp) :: xhat, yhat, zhat
347      real (kind=dp) :: dudx, dudy, dudz
348      real (kind=dp) :: scale, sc2, bigR, switcher, dswitcher
349 <
351 <    use_damped_wolf = .false.
352 <    use_undamped_wolf = .false.
353 <    if (corrMethod .eq. 1) then
354 <       use_undamped_wolf = .true.
355 <    elseif (corrMethod .eq. 2) then
356 <       use_damped_wolf = .true.
357 <    endif
349 >    real (kind=dp) :: rcuti2, rcuti3, rcuti4
350  
351      if (.not.allocated(ElectrostaticMap)) then
352         call handleError("electrostatic", "no ElectrostaticMap was present before first call of do_electrostatic_pair!")
# Line 377 | Line 369 | contains
369      yhat = d(2) * riji
370      zhat = d(3) * riji
371  
372 +    rcuti2 = rcuti*rcuti
373 +    rcuti3 = rcuti2*rcuti
374 +    rcuti4 = rcuti2*rcuti2
375 +
376 +    swi = 1.0d0 / sw
377 +
378      !! logicals
379      i_is_Charge = ElectrostaticMap(me1)%is_Charge
380      i_is_Dipole = ElectrostaticMap(me1)%is_Dipole
# Line 520 | Line 518 | contains
518  
519         if (j_is_Charge) then
520  
521 <          vterm = pre11 * q_i * q_j * riji
522 <          vpair = vpair + vterm
525 <          epot = epot + sw*vterm
521 >          if (corrMethod .eq. 1) then
522 >             vterm = pre11 * q_i * q_j * (riji - rcuti)
523  
524 <          dudr  = - sw * vterm * riji
524 >             vpair = vpair + vterm
525 >             epot = epot + sw * vterm
526 >            
527 >             dudr  = - sw * pre11 * q_i * q_j * (riji*riji*riji - rcuti2*rcuti)
528 >            
529 >             dudx = dudx + dudr * d(1)
530 >             dudy = dudy + dudr * d(2)
531 >             dudz = dudz + dudr * d(3)
532  
533 <          dudx = dudx + dudr * xhat
534 <          dudy = dudy + dudr * yhat
531 <          dudz = dudz + dudr * zhat
533 >          else
534 >             vterm = pre11 * q_i * q_j * riji
535  
536 +             vpair = vpair + vterm
537 +             epot = epot + sw * vterm
538 +            
539 +             dudr  = - sw * vterm * riji
540 +            
541 +             dudx = dudx + dudr * xhat
542 +             dudy = dudy + dudr * yhat
543 +             dudz = dudz + dudr * zhat
544 +
545 +          endif
546 +
547         endif
548  
549         if (j_is_Dipole) then
550  
551 <          if (j_is_SplitDipole) then
538 <             BigR = sqrt(r2 + 0.25_dp * d_j * d_j)
539 <             ri = 1.0_dp / BigR
540 <             scale = rij * ri
541 <          else
542 <             ri = riji
543 <             scale = 1.0_dp
544 <          endif
551 >          pref = sw * pre12 * q_i * mu_j
552  
553 <          ri2 = ri * ri
554 <          ri3 = ri2 * ri
555 <          sc2 = scale * scale
553 >          if (corrMethod .eq. 1) then
554 >             ri2 = riji * riji
555 >             ri3 = ri2 * riji
556  
557 <          pref = pre12 * q_i * mu_j
558 <          vterm = - pref * ct_j * ri2 * scale
559 <          vpair = vpair + vterm
560 <          epot = epot + sw * vterm
557 >             vterm = - pref * ct_j * (ri2 - rcuti2)
558 >             vpair = vpair + swi*vterm
559 >             epot = epot + vterm
560 >            
561 >             !! this has a + sign in the () because the rij vector is
562 >             !! r_j - r_i and the charge-dipole potential takes the origin
563 >             !! as the point dipole, which is atom j in this case.
564 >            
565 >             dudx = dudx - pref * ( ri3*( uz_j(1) - 3.0d0*ct_j*xhat) &
566 >                  - rcuti3*( uz_j(1) - 3.0d0*ct_j*d(1)*rcuti ) )
567 >             dudy = dudy - pref * ( ri3*( uz_j(2) - 3.0d0*ct_j*yhat) &
568 >                  - rcuti3*( uz_j(2) - 3.0d0*ct_j*d(2)*rcuti ) )
569 >             dudz = dudz - pref * ( ri3*( uz_j(3) - 3.0d0*ct_j*zhat) &
570 >                  - rcuti3*( uz_j(3) - 3.0d0*ct_j*d(3)*rcuti ) )
571 >            
572 >             duduz_j(1) = duduz_j(1) - pref*( ri2*xhat - d(1)*rcuti3 )
573 >             duduz_j(2) = duduz_j(2) - pref*( ri2*yhat - d(2)*rcuti3 )
574 >             duduz_j(3) = duduz_j(3) - pref*( ri2*zhat - d(3)*rcuti3 )
575  
576 <          !! this has a + sign in the () because the rij vector is
577 <          !! r_j - r_i and the charge-dipole potential takes the origin
578 <          !! as the point dipole, which is atom j in this case.
576 >          else
577 >             if (j_is_SplitDipole) then
578 >                BigR = sqrt(r2 + 0.25_dp * d_j * d_j)
579 >                ri = 1.0_dp / BigR
580 >                scale = rij * ri
581 >             else
582 >                ri = riji
583 >                scale = 1.0_dp
584 >             endif
585 >            
586 >             ri2 = ri * ri
587 >             ri3 = ri2 * ri
588 >             sc2 = scale * scale
589 >            
590 >             vterm = - pref * ct_j * ri2 * scale
591 >             vpair = vpair + swi * vterm
592 >             epot = epot + vterm
593 >            
594 >             !! this has a + sign in the () because the rij vector is
595 >             !! r_j - r_i and the charge-dipole potential takes the origin
596 >             !! as the point dipole, which is atom j in this case.
597 >            
598 >             dudx = dudx - pref * ri3 * ( uz_j(1) - 3.0d0*ct_j*xhat*sc2)
599 >             dudy = dudy - pref * ri3 * ( uz_j(2) - 3.0d0*ct_j*yhat*sc2)
600 >             dudz = dudz - pref * ri3 * ( uz_j(3) - 3.0d0*ct_j*zhat*sc2)
601 >            
602 >             duduz_j(1) = duduz_j(1) - pref * ri2 * xhat * scale
603 >             duduz_j(2) = duduz_j(2) - pref * ri2 * yhat * scale
604 >             duduz_j(3) = duduz_j(3) - pref * ri2 * zhat * scale
605  
606 <          dudx = dudx - pref * sw * ri3 * ( uz_j(1) - 3.0d0*ct_j*xhat*sc2)
560 <          dudy = dudy - pref * sw * ri3 * ( uz_j(2) - 3.0d0*ct_j*yhat*sc2)
561 <          dudz = dudz - pref * sw * ri3 * ( uz_j(3) - 3.0d0*ct_j*zhat*sc2)
562 <
563 <          duduz_j(1) = duduz_j(1) - pref * sw * ri2 * xhat * scale
564 <          duduz_j(2) = duduz_j(2) - pref * sw * ri2 * yhat * scale
565 <          duduz_j(3) = duduz_j(3) - pref * sw * ri2 * zhat * scale
566 <
606 >          endif
607         endif
608  
609         if (j_is_Quadrupole) then
# Line 575 | Line 615 | contains
615            cz2 = cz_j * cz_j
616  
617  
618 <          pref =  pre14 * q_i / 3.0_dp
579 <          vterm = pref * ri3 * (qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
580 <               qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
581 <               qzz_j * (3.0_dp*cz2 - 1.0_dp))
582 <          vpair = vpair + vterm
583 <          epot = epot + sw * vterm
618 >          pref =  sw * pre14 * q_i / 3.0_dp
619  
620 <          dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + pref * sw * ri4 * ( &
621 <               qxx_j*(6.0_dp*cx_j*ux_j(1) - 2.0_dp*xhat) + &
622 <               qyy_j*(6.0_dp*cy_j*uy_j(1) - 2.0_dp*xhat) + &
623 <               qzz_j*(6.0_dp*cz_j*uz_j(1) - 2.0_dp*xhat) )
624 <          dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + pref * sw * ri4 * ( &
625 <               qxx_j*(6.0_dp*cx_j*ux_j(2) - 2.0_dp*yhat) + &
626 <               qyy_j*(6.0_dp*cy_j*uy_j(2) - 2.0_dp*yhat) + &
627 <               qzz_j*(6.0_dp*cz_j*uz_j(2) - 2.0_dp*yhat) )
628 <          dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + pref * sw * ri4 * ( &
629 <               qxx_j*(6.0_dp*cx_j*ux_j(3) - 2.0_dp*zhat) + &
630 <               qyy_j*(6.0_dp*cy_j*uy_j(3) - 2.0_dp*zhat) + &
631 <               qzz_j*(6.0_dp*cz_j*uz_j(3) - 2.0_dp*zhat) )
632 <
633 <          dudux_j(1) = dudux_j(1) + pref * sw * ri3 * (qxx_j*6.0_dp*cx_j*xhat)
634 <          dudux_j(2) = dudux_j(2) + pref * sw * ri3 * (qxx_j*6.0_dp*cx_j*yhat)
635 <          dudux_j(3) = dudux_j(3) + pref * sw * ri3 * (qxx_j*6.0_dp*cx_j*zhat)
636 <
637 <          duduy_j(1) = duduy_j(1) + pref * sw * ri3 * (qyy_j*6.0_dp*cy_j*xhat)
638 <          duduy_j(2) = duduy_j(2) + pref * sw * ri3 * (qyy_j*6.0_dp*cy_j*yhat)
639 <          duduy_j(3) = duduy_j(3) + pref * sw * ri3 * (qyy_j*6.0_dp*cy_j*zhat)
640 <
641 <          duduz_j(1) = duduz_j(1) + pref * sw * ri3 * (qzz_j*6.0_dp*cz_j*xhat)
642 <          duduz_j(2) = duduz_j(2) + pref * sw * ri3 * (qzz_j*6.0_dp*cz_j*yhat)
643 <          duduz_j(3) = duduz_j(3) + pref * sw * ri3 * (qzz_j*6.0_dp*cz_j*zhat)
620 >          if (corrMethod .eq. 1) then
621 >             vterm1 = pref * ri3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
622 >                  qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
623 >                  qzz_j * (3.0_dp*cz2 - 1.0_dp) )
624 >             vterm2 = pref * rcuti3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
625 >                  qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
626 >                  qzz_j * (3.0_dp*cz2 - 1.0_dp) )
627 >             vpair = vpair + swi*( vterm1 - vterm2 )
628 >             epot = epot + ( vterm1 - vterm2 )
629 >            
630 >             dudx = dudx - (5.0_dp * &
631 >                  (vterm1*riji*xhat - vterm2*rcuti2*d(1))) + pref * ( &
632 >                  (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(1)) - &
633 >                  qxx_j*2.0_dp*(xhat - rcuti*d(1))) + &
634 >                  (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(1)) - &
635 >                  qyy_j*2.0_dp*(xhat - rcuti*d(1))) + &
636 >                  (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(1)) - &
637 >                  qzz_j*2.0_dp*(xhat - rcuti*d(1))) )
638 >             dudy = dudy - (5.0_dp * &
639 >                  (vterm1*riji*yhat - vterm2*rcuti2*d(2))) + pref * ( &
640 >                  (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(2)) - &
641 >                  qxx_j*2.0_dp*(yhat - rcuti*d(2))) + &
642 >                  (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(2)) - &
643 >                  qyy_j*2.0_dp*(yhat - rcuti*d(2))) + &
644 >                  (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(2)) - &
645 >                  qzz_j*2.0_dp*(yhat - rcuti*d(2))) )
646 >             dudz = dudz - (5.0_dp * &
647 >                  (vterm1*riji*zhat - vterm2*rcuti2*d(3))) + pref * ( &
648 >                  (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(3)) - &
649 >                  qxx_j*2.0_dp*(zhat - rcuti*d(3))) + &
650 >                  (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(3)) - &
651 >                  qyy_j*2.0_dp*(zhat - rcuti*d(3))) + &
652 >                  (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(3)) - &
653 >                  qzz_j*2.0_dp*(zhat - rcuti*d(3))) )
654 >            
655 >             dudux_j(1) = dudux_j(1) + pref * (ri3*(qxx_j*6.0_dp*cx_j*xhat) - &
656 >                  rcuti4*(qxx_j*6.0_dp*cx_j*d(1)))
657 >             dudux_j(2) = dudux_j(2) + pref * (ri3*(qxx_j*6.0_dp*cx_j*yhat) - &
658 >                  rcuti4*(qxx_j*6.0_dp*cx_j*d(2)))
659 >             dudux_j(3) = dudux_j(3) + pref * (ri3*(qxx_j*6.0_dp*cx_j*zhat) - &
660 >                  rcuti4*(qxx_j*6.0_dp*cx_j*d(3)))
661 >            
662 >             duduy_j(1) = duduy_j(1) + pref * (ri3*(qyy_j*6.0_dp*cy_j*xhat) - &
663 >                  rcuti4*(qyy_j*6.0_dp*cx_j*d(1)))
664 >             duduy_j(2) = duduy_j(2) + pref * (ri3*(qyy_j*6.0_dp*cy_j*yhat) - &
665 >                  rcuti4*(qyy_j*6.0_dp*cx_j*d(2)))
666 >             duduy_j(3) = duduy_j(3) + pref * (ri3*(qyy_j*6.0_dp*cy_j*zhat) - &
667 >                  rcuti4*(qyy_j*6.0_dp*cx_j*d(3)))
668 >            
669 >             duduz_j(1) = duduz_j(1) + pref * (ri3*(qzz_j*6.0_dp*cz_j*xhat) - &
670 >                  rcuti4*(qzz_j*6.0_dp*cx_j*d(1)))
671 >             duduz_j(2) = duduz_j(2) + pref * (ri3*(qzz_j*6.0_dp*cz_j*yhat) - &
672 >                  rcuti4*(qzz_j*6.0_dp*cx_j*d(2)))
673 >             duduz_j(3) = duduz_j(3) + pref * (ri3*(qzz_j*6.0_dp*cz_j*zhat) - &
674 >                  rcuti4*(qzz_j*6.0_dp*cx_j*d(3)))
675 >        
676 >          else
677 >             vterm = pref * ri3 * (qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
678 >                  qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
679 >                  qzz_j * (3.0_dp*cz2 - 1.0_dp))
680 >             vpair = vpair + swi * vterm
681 >             epot = epot + vterm
682 >            
683 >             dudx = dudx - 5.0_dp*vterm*riji*xhat + pref * ri4 * ( &
684 >                  qxx_j*(6.0_dp*cx_j*ux_j(1) - 2.0_dp*xhat) + &
685 >                  qyy_j*(6.0_dp*cy_j*uy_j(1) - 2.0_dp*xhat) + &
686 >                  qzz_j*(6.0_dp*cz_j*uz_j(1) - 2.0_dp*xhat) )
687 >             dudy = dudy - 5.0_dp*vterm*riji*yhat + pref * ri4 * ( &
688 >                  qxx_j*(6.0_dp*cx_j*ux_j(2) - 2.0_dp*yhat) + &
689 >                  qyy_j*(6.0_dp*cy_j*uy_j(2) - 2.0_dp*yhat) + &
690 >                  qzz_j*(6.0_dp*cz_j*uz_j(2) - 2.0_dp*yhat) )
691 >             dudz = dudz - 5.0_dp*vterm*riji*zhat + pref * ri4 * ( &
692 >                  qxx_j*(6.0_dp*cx_j*ux_j(3) - 2.0_dp*zhat) + &
693 >                  qyy_j*(6.0_dp*cy_j*uy_j(3) - 2.0_dp*zhat) + &
694 >                  qzz_j*(6.0_dp*cz_j*uz_j(3) - 2.0_dp*zhat) )
695 >            
696 >             dudux_j(1) = dudux_j(1) + pref * ri3*(qxx_j*6.0_dp*cx_j*xhat)
697 >             dudux_j(2) = dudux_j(2) + pref * ri3*(qxx_j*6.0_dp*cx_j*yhat)
698 >             dudux_j(3) = dudux_j(3) + pref * ri3*(qxx_j*6.0_dp*cx_j*zhat)
699 >            
700 >             duduy_j(1) = duduy_j(1) + pref * ri3*(qyy_j*6.0_dp*cy_j*xhat)
701 >             duduy_j(2) = duduy_j(2) + pref * ri3*(qyy_j*6.0_dp*cy_j*yhat)
702 >             duduy_j(3) = duduy_j(3) + pref * ri3*(qyy_j*6.0_dp*cy_j*zhat)
703 >            
704 >             duduz_j(1) = duduz_j(1) + pref * ri3*(qzz_j*6.0_dp*cz_j*xhat)
705 >             duduz_j(2) = duduz_j(2) + pref * ri3*(qzz_j*6.0_dp*cz_j*yhat)
706 >             duduz_j(3) = duduz_j(3) + pref * ri3*(qzz_j*6.0_dp*cz_j*zhat)
707 >          
708 >          endif
709         endif
610 –
710      endif
711  
712      if (i_is_Dipole) then
713  
714         if (j_is_Charge) then
715  
716 <          if (i_is_SplitDipole) then
618 <             BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
619 <             ri = 1.0_dp / BigR
620 <             scale = rij * ri
621 <          else
622 <             ri = riji
623 <             scale = 1.0_dp
624 <          endif
716 >          pref = sw * pre12 * q_j * mu_i
717  
718 <          ri2 = ri * ri
719 <          ri3 = ri2 * ri
720 <          sc2 = scale * scale
718 >          if (corrMethod .eq. 1) then
719 >             ri2 = riji * riji
720 >             ri3 = ri2 * riji
721  
722 <          pref = pre12 * q_j * mu_i
723 <          vterm = pref * ct_i * ri2 * scale
724 <          vpair = vpair + vterm
725 <          epot = epot + sw * vterm
726 <
727 <          dudx = dudx + pref * sw * ri3 * ( uz_i(1) - 3.0d0 * ct_i * xhat*sc2)
728 <          dudy = dudy + pref * sw * ri3 * ( uz_i(2) - 3.0d0 * ct_i * yhat*sc2)
729 <          dudz = dudz + pref * sw * ri3 * ( uz_i(3) - 3.0d0 * ct_i * zhat*sc2)
730 <
731 <          duduz_i(1) = duduz_i(1) + pref * sw * ri2 * xhat * scale
732 <          duduz_i(2) = duduz_i(2) + pref * sw * ri2 * yhat * scale
733 <          duduz_i(3) = duduz_i(3) + pref * sw * ri2 * zhat * scale
734 <       endif
735 <
736 <       if (j_is_Dipole) then
722 >             vterm = pref * ct_i * (ri2 - rcuti2)
723 >             vpair = vpair + swi * vterm
724 >             epot = epot + vterm
725 >            
726 >             !! this has a + sign in the () because the rij vector is
727 >             !! r_j - r_i and the charge-dipole potential takes the origin
728 >             !! as the point dipole, which is atom j in this case.
729 >            
730 >             dudx = dudx + pref * ( ri3*( uz_i(1) - 3.0d0*ct_i*xhat) &
731 >                  - rcuti3*( uz_i(1) - 3.0d0*ct_i*d(1)*rcuti ) )
732 >             dudy = dudy + pref * ( ri3*( uz_i(2) - 3.0d0*ct_i*yhat) &
733 >                  - rcuti3*( uz_i(2) - 3.0d0*ct_i*d(2)*rcuti ) )
734 >             dudz = dudz + pref * ( ri3*( uz_i(3) - 3.0d0*ct_i*zhat) &
735 >                  - rcuti3*( uz_i(3) - 3.0d0*ct_i*d(3)*rcuti ) )
736 >            
737 >             duduz_i(1) = duduz_i(1) - pref*( ri2*xhat - d(1)*rcuti3 )
738 >             duduz_i(2) = duduz_i(2) - pref*( ri2*yhat - d(2)*rcuti3 )
739 >             duduz_i(3) = duduz_i(3) - pref*( ri2*zhat - d(3)*rcuti3 )
740  
646 –          if (i_is_SplitDipole) then
647 –             if (j_is_SplitDipole) then
648 –                BigR = sqrt(r2 + 0.25_dp * d_i * d_i + 0.25_dp * d_j * d_j)
649 –             else
650 –                BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
651 –             endif
652 –             ri = 1.0_dp / BigR
653 –             scale = rij * ri                
741            else
742 <             if (j_is_SplitDipole) then
743 <                BigR = sqrt(r2 + 0.25_dp * d_j * d_j)
742 >             if (i_is_SplitDipole) then
743 >                BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
744                  ri = 1.0_dp / BigR
745 <                scale = rij * ri                            
746 <             else                
745 >                scale = rij * ri
746 >             else
747                  ri = riji
748                  scale = 1.0_dp
749               endif
750 +            
751 +             ri2 = ri * ri
752 +             ri3 = ri2 * ri
753 +             sc2 = scale * scale
754 +            
755 +             vterm = pref * ct_i * ri2 * scale
756 +             vpair = vpair + swi * vterm
757 +             epot = epot + vterm
758 +            
759 +             dudx = dudx + pref * ri3 * ( uz_i(1) - 3.0d0 * ct_i * xhat*sc2)
760 +             dudy = dudy + pref * ri3 * ( uz_i(2) - 3.0d0 * ct_i * yhat*sc2)
761 +             dudz = dudz + pref * ri3 * ( uz_i(3) - 3.0d0 * ct_i * zhat*sc2)
762 +            
763 +             duduz_i(1) = duduz_i(1) + pref * ri2 * xhat * scale
764 +             duduz_i(2) = duduz_i(2) + pref * ri2 * yhat * scale
765 +             duduz_i(3) = duduz_i(3) + pref * ri2 * zhat * scale
766            endif
767 +       endif
768  
769 <          ct_ij = uz_i(1)*uz_j(1) + uz_i(2)*uz_j(2) + uz_i(3)*uz_j(3)
769 >       if (j_is_Dipole) then
770  
771 <          ri2 = ri * ri
668 <          ri3 = ri2 * ri
669 <          ri4 = ri2 * ri2
670 <          sc2 = scale * scale
771 >          pref = sw * pre22 * mu_i * mu_j
772  
773 <          pref = pre22 * mu_i * mu_j
774 <          vterm = pref * ri3 * (ct_ij - 3.0d0 * ct_i * ct_j * sc2)
775 <          vpair = vpair + vterm
776 <          epot = epot + sw * vterm
773 >          if (corrMethod .eq. 1) then
774 >             ri2 = riji * riji
775 >             ri3 = ri2 * riji
776 >             ri4 = ri2 * ri2
777  
778 <          a1 = 5.0d0 * ct_i * ct_j * sc2 - ct_ij
779 <
780 <          dudx = dudx + pref*sw*3.0d0*ri4*scale &
781 <                         *(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
782 <          dudy = dudy + pref*sw*3.0d0*ri4*scale &
783 <                         *(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
784 <          dudz = dudz + pref*sw*3.0d0*ri4*scale &
785 <                         *(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
786 <
787 <          duduz_i(1) = duduz_i(1) + pref*sw*ri3 &
788 <                                     *(uz_j(1) - 3.0d0*ct_j*xhat*sc2)
789 <          duduz_i(2) = duduz_i(2) + pref*sw*ri3 &
790 <                                     *(uz_j(2) - 3.0d0*ct_j*yhat*sc2)
791 <          duduz_i(3) = duduz_i(3) + pref*sw*ri3 &
792 <                                     *(uz_j(3) - 3.0d0*ct_j*zhat*sc2)
793 <
794 <          duduz_j(1) = duduz_j(1) + pref*sw*ri3 &
795 <                                     *(uz_i(1) - 3.0d0*ct_i*xhat*sc2)
796 <          duduz_j(2) = duduz_j(2) + pref*sw*ri3 &
797 <                                     *(uz_i(2) - 3.0d0*ct_i*yhat*sc2)
798 <          duduz_j(3) = duduz_j(3) + pref*sw*ri3 &
799 <                                     *(uz_i(3) - 3.0d0*ct_i*zhat*sc2)
778 >             vterm = pref * (ri3 - rcuti3) * (ct_ij - 3.0d0 * ct_i * ct_j)
779 >             vpair = vpair + swi * vterm
780 >             epot = epot + vterm
781 >            
782 >             a1 = 5.0d0 * ct_i * ct_j - ct_ij
783 >            
784 >             dudx = dudx + pref*3.0d0*ri4 &
785 >                  *(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1)) - &
786 >                  pref*3.0d0*rcuti4*(a1*rcuti*d(1)-ct_i*uz_j(1)-ct_j*uz_i(1))
787 >             dudy = dudy + pref*3.0d0*ri4 &
788 >                  *(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2)) - &
789 >                  pref*3.0d0*rcuti4*(a1*rcuti*d(2)-ct_i*uz_j(2)-ct_j*uz_i(2))
790 >             dudz = dudz + pref*3.0d0*ri4 &
791 >                  *(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3)) - &
792 >                  pref*3.0d0*rcuti4*(a1*rcuti*d(3)-ct_i*uz_j(3)-ct_j*uz_i(3))
793 >            
794 >             duduz_i(1) = duduz_i(1) + pref*(ri3*(uz_j(1) - 3.0d0*ct_j*xhat) &
795 >                  - rcuti3*(uz_j(1) - 3.0d0*ct_j*d(1)*rcuti))
796 >             duduz_i(2) = duduz_i(2) + pref*(ri3*(uz_j(2) - 3.0d0*ct_j*yhat) &
797 >                  - rcuti3*(uz_j(2) - 3.0d0*ct_j*d(2)*rcuti))
798 >             duduz_i(3) = duduz_i(3) + pref*(ri3*(uz_j(3) - 3.0d0*ct_j*zhat) &
799 >                  - rcuti3*(uz_j(3) - 3.0d0*ct_j*d(3)*rcuti))
800 >             duduz_j(1) = duduz_j(1) + pref*(ri3*(uz_i(1) - 3.0d0*ct_i*xhat) &
801 >                  - rcuti3*(uz_i(1) - 3.0d0*ct_i*d(1)*rcuti))
802 >             duduz_j(2) = duduz_j(2) + pref*(ri3*(uz_i(2) - 3.0d0*ct_i*yhat) &
803 >                  - rcuti3*(uz_i(2) - 3.0d0*ct_i*d(2)*rcuti))
804 >             duduz_j(3) = duduz_j(3) + pref*(ri3*(uz_i(3) - 3.0d0*ct_i*zhat) &
805 >                  - rcuti3*(uz_i(3) - 3.0d0*ct_i*d(3)*rcuti))
806 >          else
807 >            
808 >             if (i_is_SplitDipole) then
809 >                if (j_is_SplitDipole) then
810 >                   BigR = sqrt(r2 + 0.25_dp * d_i * d_i + 0.25_dp * d_j * d_j)
811 >                else
812 >                   BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
813 >                endif
814 >                ri = 1.0_dp / BigR
815 >                scale = rij * ri                
816 >             else
817 >                if (j_is_SplitDipole) then
818 >                   BigR = sqrt(r2 + 0.25_dp * d_j * d_j)
819 >                   ri = 1.0_dp / BigR
820 >                   scale = rij * ri                            
821 >                else                
822 >                   ri = riji
823 >                   scale = 1.0_dp
824 >                endif
825 >             endif
826 >            
827 >             ct_ij = uz_i(1)*uz_j(1) + uz_i(2)*uz_j(2) + uz_i(3)*uz_j(3)
828 >            
829 >             ri2 = ri * ri
830 >             ri3 = ri2 * ri
831 >             ri4 = ri2 * ri2
832 >             sc2 = scale * scale
833 >            
834 >             vterm = pref * ri3 * (ct_ij - 3.0d0 * ct_i * ct_j * sc2)
835 >             vpair = vpair + swi * vterm
836 >             epot = epot + vterm
837 >            
838 >             a1 = 5.0d0 * ct_i * ct_j * sc2 - ct_ij
839 >            
840 >             dudx = dudx + pref*3.0d0*ri4*scale &
841 >                  *(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
842 >             dudy = dudy + pref*3.0d0*ri4*scale &
843 >                  *(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
844 >             dudz = dudz + pref*3.0d0*ri4*scale &
845 >                  *(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
846 >            
847 >             duduz_i(1) = duduz_i(1) + pref*ri3 &
848 >                  *(uz_j(1) - 3.0d0*ct_j*xhat*sc2)
849 >             duduz_i(2) = duduz_i(2) + pref*ri3 &
850 >                  *(uz_j(2) - 3.0d0*ct_j*yhat*sc2)
851 >             duduz_i(3) = duduz_i(3) + pref*ri3 &
852 >                  *(uz_j(3) - 3.0d0*ct_j*zhat*sc2)
853 >            
854 >             duduz_j(1) = duduz_j(1) + pref*ri3 &
855 >                  *(uz_i(1) - 3.0d0*ct_i*xhat*sc2)
856 >             duduz_j(2) = duduz_j(2) + pref*ri3 &
857 >                  *(uz_i(2) - 3.0d0*ct_i*yhat*sc2)
858 >             duduz_j(3) = duduz_j(3) + pref*ri3 &
859 >                  *(uz_i(3) - 3.0d0*ct_i*zhat*sc2)
860 >          endif
861         endif
700 –
862      endif
863  
864      if (i_is_Quadrupole) then
# Line 710 | Line 871 | contains
871            cy2 = cy_i * cy_i
872            cz2 = cz_i * cz_i
873  
874 <          pref = pre14 * q_j / 3.0_dp
714 <          vterm = pref * ri3 * (qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
715 <               qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
716 <               qzz_i * (3.0_dp*cz2 - 1.0_dp))
717 <          vpair = vpair + vterm
718 <          epot = epot + sw * vterm
874 >          pref = sw * pre14 * q_j / 3.0_dp
875  
876 <          dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + pref * sw * ri4 * ( &
877 <               qxx_i*(6.0_dp*cx_i*ux_i(1) - 2.0_dp*xhat) + &
878 <               qyy_i*(6.0_dp*cy_i*uy_i(1) - 2.0_dp*xhat) + &
879 <               qzz_i*(6.0_dp*cz_i*uz_i(1) - 2.0_dp*xhat) )
880 <          dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + pref * sw * ri4 * ( &
881 <               qxx_i*(6.0_dp*cx_i*ux_i(2) - 2.0_dp*yhat) + &
882 <               qyy_i*(6.0_dp*cy_i*uy_i(2) - 2.0_dp*yhat) + &
883 <               qzz_i*(6.0_dp*cz_i*uz_i(2) - 2.0_dp*yhat) )
884 <          dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + pref * sw * ri4 * ( &
885 <               qxx_i*(6.0_dp*cx_i*ux_i(3) - 2.0_dp*zhat) + &
886 <               qyy_i*(6.0_dp*cy_i*uy_i(3) - 2.0_dp*zhat) + &
887 <               qzz_i*(6.0_dp*cz_i*uz_i(3) - 2.0_dp*zhat) )
876 >          if (corrMethod .eq. 1) then
877 >             vterm1 = pref * ri3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
878 >                  qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
879 >                  qzz_i * (3.0_dp*cz2 - 1.0_dp) )
880 >             vterm2 = pref * rcuti3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
881 >                  qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
882 >                  qzz_i * (3.0_dp*cz2 - 1.0_dp) )
883 >             vpair = vpair + swi * ( vterm1 - vterm2 )
884 >             epot = epot + ( vterm1 - vterm2 )
885 >            
886 >             dudx = dudx - (5.0_dp*(vterm1*riji*xhat - vterm2*rcuti2*d(1))) + &
887 >                  pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(1)) - &
888 >                  qxx_i*2.0_dp*(xhat - rcuti*d(1))) + &
889 >                  (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(1)) - &
890 >                  qyy_i*2.0_dp*(xhat - rcuti*d(1))) + &
891 >                  (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(1)) - &
892 >                  qzz_i*2.0_dp*(xhat - rcuti*d(1))) )
893 >             dudy = dudy - (5.0_dp*(vterm1*riji*yhat - vterm2*rcuti2*d(2))) + &
894 >                  pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(2)) - &
895 >                  qxx_i*2.0_dp*(yhat - rcuti*d(2))) + &
896 >                  (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(2)) - &
897 >                  qyy_i*2.0_dp*(yhat - rcuti*d(2))) + &
898 >                  (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(2)) - &
899 >                  qzz_i*2.0_dp*(yhat - rcuti*d(2))) )
900 >             dudz = dudz - (5.0_dp*(vterm1*riji*zhat - vterm2*rcuti2*d(3))) + &
901 >                  pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(3)) - &
902 >                  qxx_i*2.0_dp*(zhat - rcuti*d(3))) + &
903 >                  (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(3)) - &
904 >                  qyy_i*2.0_dp*(zhat - rcuti*d(3))) + &
905 >                  (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(3)) - &
906 >                  qzz_i*2.0_dp*(zhat - rcuti*d(3))) )
907 >            
908 >             dudux_i(1) = dudux_i(1) + pref * (ri3*(qxx_i*6.0_dp*cx_i*xhat) - &
909 >                  rcuti4*(qxx_i*6.0_dp*cx_i*d(1)))
910 >             dudux_i(2) = dudux_i(2) + pref * (ri3*(qxx_i*6.0_dp*cx_i*yhat) - &
911 >                  rcuti4*(qxx_i*6.0_dp*cx_i*d(2)))
912 >             dudux_i(3) = dudux_i(3) + pref * (ri3*(qxx_i*6.0_dp*cx_i*zhat) - &
913 >                  rcuti4*(qxx_i*6.0_dp*cx_i*d(3)))
914 >            
915 >             duduy_i(1) = duduy_i(1) + pref * (ri3*(qyy_i*6.0_dp*cy_i*xhat) - &
916 >                  rcuti4*(qyy_i*6.0_dp*cx_i*d(1)))
917 >             duduy_i(2) = duduy_i(2) + pref * (ri3*(qyy_i*6.0_dp*cy_i*yhat) - &
918 >                  rcuti4*(qyy_i*6.0_dp*cx_i*d(2)))
919 >             duduy_i(3) = duduy_i(3) + pref * (ri3*(qyy_i*6.0_dp*cy_i*zhat) - &
920 >                  rcuti4*(qyy_i*6.0_dp*cx_i*d(3)))
921 >            
922 >             duduz_i(1) = duduz_i(1) + pref * (ri3*(qzz_i*6.0_dp*cz_i*xhat) - &
923 >                  rcuti4*(qzz_i*6.0_dp*cx_i*d(1)))
924 >             duduz_i(2) = duduz_i(2) + pref * (ri3*(qzz_i*6.0_dp*cz_i*yhat) - &
925 >                  rcuti4*(qzz_i*6.0_dp*cx_i*d(2)))
926 >             duduz_i(3) = duduz_i(3) + pref * (ri3*(qzz_i*6.0_dp*cz_i*zhat) - &
927 >                  rcuti4*(qzz_i*6.0_dp*cx_i*d(3)))
928  
929 <          dudux_i(1) = dudux_i(1) + pref * sw * ri3 * (qxx_i*6.0_dp*cx_i*xhat)
930 <          dudux_i(2) = dudux_i(2) + pref * sw * ri3 * (qxx_i*6.0_dp*cx_i*yhat)
931 <          dudux_i(3) = dudux_i(3) + pref * sw * ri3 * (qxx_i*6.0_dp*cx_i*zhat)
932 <
933 <          duduy_i(1) = duduy_i(1) + pref * sw * ri3 * (qyy_i*6.0_dp*cy_i*xhat)
934 <          duduy_i(2) = duduy_i(2) + pref * sw * ri3 * (qyy_i*6.0_dp*cy_i*yhat)
935 <          duduy_i(3) = duduy_i(3) + pref * sw * ri3 * (qyy_i*6.0_dp*cy_i*zhat)
936 <
937 <          duduz_i(1) = duduz_i(1) + pref * sw * ri3 * (qzz_i*6.0_dp*cz_i*xhat)
938 <          duduz_i(2) = duduz_i(2) + pref * sw * ri3 * (qzz_i*6.0_dp*cz_i*yhat)
939 <          duduz_i(3) = duduz_i(3) + pref * sw * ri3 * (qzz_i*6.0_dp*cz_i*zhat)
929 >          else
930 >             vterm = pref * ri3 * (qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
931 >                  qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
932 >                  qzz_i * (3.0_dp*cz2 - 1.0_dp))
933 >             vpair = vpair + swi * vterm
934 >             epot = epot + vterm
935 >            
936 >             dudx = dudx - 5.0_dp*vterm*riji*xhat + pref * ri4 * ( &
937 >                  qxx_i*(6.0_dp*cx_i*ux_i(1) - 2.0_dp*xhat) + &
938 >                  qyy_i*(6.0_dp*cy_i*uy_i(1) - 2.0_dp*xhat) + &
939 >                  qzz_i*(6.0_dp*cz_i*uz_i(1) - 2.0_dp*xhat) )
940 >             dudy = dudy - 5.0_dp*vterm*riji*yhat + pref * ri4 * ( &
941 >                  qxx_i*(6.0_dp*cx_i*ux_i(2) - 2.0_dp*yhat) + &
942 >                  qyy_i*(6.0_dp*cy_i*uy_i(2) - 2.0_dp*yhat) + &
943 >                  qzz_i*(6.0_dp*cz_i*uz_i(2) - 2.0_dp*yhat) )
944 >             dudz = dudz - 5.0_dp*vterm*riji*zhat + pref * ri4 * ( &
945 >                  qxx_i*(6.0_dp*cx_i*ux_i(3) - 2.0_dp*zhat) + &
946 >                  qyy_i*(6.0_dp*cy_i*uy_i(3) - 2.0_dp*zhat) + &
947 >                  qzz_i*(6.0_dp*cz_i*uz_i(3) - 2.0_dp*zhat) )
948 >            
949 >             dudux_i(1) = dudux_i(1) + pref * ri3*(qxx_i*6.0_dp*cx_i*xhat)
950 >             dudux_i(2) = dudux_i(2) + pref * ri3*(qxx_i*6.0_dp*cx_i*yhat)
951 >             dudux_i(3) = dudux_i(3) + pref * ri3*(qxx_i*6.0_dp*cx_i*zhat)
952 >            
953 >             duduy_i(1) = duduy_i(1) + pref * ri3*(qyy_i*6.0_dp*cy_i*xhat)
954 >             duduy_i(2) = duduy_i(2) + pref * ri3*(qyy_i*6.0_dp*cy_i*yhat)
955 >             duduy_i(3) = duduy_i(3) + pref * ri3*(qyy_i*6.0_dp*cy_i*zhat)
956 >            
957 >             duduz_i(1) = duduz_i(1) + pref * ri3*(qzz_i*6.0_dp*cz_i*xhat)
958 >             duduz_i(2) = duduz_i(2) + pref * ri3*(qzz_i*6.0_dp*cz_i*yhat)
959 >             duduz_i(3) = duduz_i(3) + pref * ri3*(qzz_i*6.0_dp*cz_i*zhat)
960 >          endif
961         endif
962      endif
963  

Diff Legend

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