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 2204 by gezelter, Fri Apr 15 22:04:00 2005 UTC vs.
Revision 2296 by chrisfen, Thu Sep 15 00:13:56 2005 UTC

# Line 85 | Line 85 | module electrostatic_module
85       logical :: is_Dipole = .false.
86       logical :: is_SplitDipole = .false.
87       logical :: is_Quadrupole = .false.
88 +     logical :: is_Tap = .false.
89       real(kind=DP) :: charge = 0.0_DP
90       real(kind=DP) :: dipole_moment = 0.0_DP
91       real(kind=DP) :: split_dipole_distance = 0.0_DP
# Line 96 | Line 97 | contains
97   contains
98  
99    subroutine newElectrostaticType(c_ident, is_Charge, is_Dipole, &
100 <       is_SplitDipole, is_Quadrupole, status)
100 >       is_SplitDipole, is_Quadrupole, is_Tap, status)
101  
102      integer, intent(in) :: c_ident
103      logical, intent(in) :: is_Charge
104      logical, intent(in) :: is_Dipole
105      logical, intent(in) :: is_SplitDipole
106      logical, intent(in) :: is_Quadrupole
107 +    logical, intent(in) :: is_Tap
108      integer, intent(out) :: status
109      integer :: nAtypes, myATID, i, j
110  
# Line 139 | Line 141 | contains
141      ElectrostaticMap(myATID)%is_Dipole = is_Dipole
142      ElectrostaticMap(myATID)%is_SplitDipole = is_SplitDipole
143      ElectrostaticMap(myATID)%is_Quadrupole = is_Quadrupole
144 +    ElectrostaticMap(myATID)%is_Tap = is_Tap
145  
146    end subroutine newElectrostaticType
147  
# Line 304 | 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)
310 >       vpair, fpair, pot, eFrame, f, t, do_pot, corrMethod, rcuti)
311  
312      logical, intent(in) :: do_pot
313  
314      integer, intent(in) :: atom1, atom2
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 328 | Line 332 | contains
332  
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      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 337 | 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
348 >    real (kind=dp) :: scale, sc2, bigR, switcher, dswitcher
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 363 | Line 369 | contains
369      yhat = d(2) * riji
370      zhat = d(3) * riji
371  
372 <    !! logicals
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
381      i_is_SplitDipole = ElectrostaticMap(me1)%is_SplitDipole
382      i_is_Quadrupole = ElectrostaticMap(me1)%is_Quadrupole
383 +    i_is_Tap = ElectrostaticMap(me1)%is_Tap
384  
385      j_is_Charge = ElectrostaticMap(me2)%is_Charge
386      j_is_Dipole = ElectrostaticMap(me2)%is_Dipole
387      j_is_SplitDipole = ElectrostaticMap(me2)%is_SplitDipole
388      j_is_Quadrupole = ElectrostaticMap(me2)%is_Quadrupole
389 +    j_is_Tap = ElectrostaticMap(me2)%is_Tap
390  
391      if (i_is_Charge) then
392         q_i = ElectrostaticMap(me1)%charge      
# Line 428 | Line 441 | contains
441         cz_i = uz_i(1)*xhat + uz_i(2)*yhat + uz_i(3)*zhat
442      endif
443  
431 –
444      if (j_is_Charge) then
445         q_j = ElectrostaticMap(me2)%charge      
446      endif
# Line 480 | Line 492 | contains
492         cy_j = uy_j(1)*xhat + uy_j(2)*yhat + uy_j(3)*zhat
493         cz_j = uz_j(1)*xhat + uz_j(2)*yhat + uz_j(3)*zhat
494      endif
495 +  
496 + !!$    switcher = 1.0d0
497 + !!$    dswitcher = 0.0d0
498 + !!$    ebalance = 0.0d0
499 + !!$    ! weaken the dipole interaction at close range for TAP water
500 + !!$    if (j_is_Tap .and. i_is_Tap) then
501 + !!$      call calc_switch(rij, mu_i, switcher, dswitcher)
502 + !!$    endif
503  
504      epot = 0.0_dp
505      dudx = 0.0_dp
# Line 498 | Line 518 | contains
518  
519         if (j_is_Charge) then
520  
521 <          vterm = pre11 * q_i * q_j * riji
522 <          vpair = vpair + vterm
503 <          epot = epot + sw*vterm
504 <
505 <          dudr  = - sw * vterm * riji
521 >          if (corrMethod .eq. 1) then
522 >             vterm = pre11 * q_i * q_j * (riji - rcuti)
523  
524 <          dudx = dudx + dudr * xhat
525 <          dudy = dudy + dudr * yhat
526 <          dudz = dudz + dudr * zhat
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 +          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
516 <             BigR = sqrt(r2 + 0.25_dp * d_j * d_j)
517 <             ri = 1.0_dp / BigR
518 <             scale = rij * ri
519 <          else
520 <             ri = riji
521 <             scale = 1.0_dp
522 <          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)
538 <          dudy = dudy - pref * sw * ri3 * ( uz_j(2) - 3.0d0*ct_j*yhat*sc2)
539 <          dudz = dudz - pref * sw * ri3 * ( uz_j(3) - 3.0d0*ct_j*zhat*sc2)
540 <
541 <          duduz_j(1) = duduz_j(1) - pref * sw * ri2 * xhat * scale
542 <          duduz_j(2) = duduz_j(2) - pref * sw * ri2 * yhat * scale
543 <          duduz_j(3) = duduz_j(3) - pref * sw * ri2 * zhat * scale
544 <
606 >          endif
607         endif
608  
609         if (j_is_Quadrupole) then
# Line 553 | Line 615 | contains
615            cz2 = cz_j * cz_j
616  
617  
618 <          pref =  pre14 * q_i / 3.0_dp
557 <          vterm = pref * ri3 * (qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
558 <               qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
559 <               qzz_j * (3.0_dp*cz2 - 1.0_dp))
560 <          vpair = vpair + vterm
561 <          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
588 –
710      endif
711  
712      if (i_is_Dipole) then
713  
714         if (j_is_Charge) then
715  
716 <          if (i_is_SplitDipole) then
596 <             BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
597 <             ri = 1.0_dp / BigR
598 <             scale = rij * ri
599 <          else
600 <             ri = riji
601 <             scale = 1.0_dp
602 <          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  
624 –          if (i_is_SplitDipole) then
625 –             if (j_is_SplitDipole) then
626 –                BigR = sqrt(r2 + 0.25_dp * d_i * d_i + 0.25_dp * d_j * d_j)
627 –             else
628 –                BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
629 –             endif
630 –             ri = 1.0_dp / BigR
631 –             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
646 <          ri3 = ri2 * ri
647 <          ri4 = ri2 * ri2
648 <          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*(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
781 <          dudy=dudy+pref*sw*3.0d0*ri4*scale*(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
782 <          dudz=dudz+pref*sw*3.0d0*ri4*scale*(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
783 <
784 <          duduz_i(1) = duduz_i(1) + pref*sw*ri3*(uz_j(1) - 3.0d0*ct_j*xhat*sc2)
785 <          duduz_i(2) = duduz_i(2) + pref*sw*ri3*(uz_j(2) - 3.0d0*ct_j*yhat*sc2)
786 <          duduz_i(3) = duduz_i(3) + pref*sw*ri3*(uz_j(3) - 3.0d0*ct_j*zhat*sc2)
787 <
788 <          duduz_j(1) = duduz_j(1) + pref*sw*ri3*(uz_i(1) - 3.0d0*ct_i*xhat*sc2)
789 <          duduz_j(2) = duduz_j(2) + pref*sw*ri3*(uz_i(2) - 3.0d0*ct_i*yhat*sc2)
790 <          duduz_j(3) = duduz_j(3) + pref*sw*ri3*(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
669 –
862      endif
863  
864      if (i_is_Quadrupole) then
# Line 679 | Line 871 | contains
871            cy2 = cy_i * cy_i
872            cz2 = cz_i * cz_i
873  
874 <          pref = pre14 * q_j / 3.0_dp
683 <          vterm = pref * ri3 * (qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
684 <               qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
685 <               qzz_i * (3.0_dp*cz2 - 1.0_dp))
686 <          vpair = vpair + vterm
687 <          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  
# Line 822 | Line 1070 | contains
1070      return
1071    end subroutine doElectrostaticPair
1072  
1073 +  !! calculates the switching functions and their derivatives for a given
1074 +  subroutine calc_switch(r, mu, scale, dscale)
1075  
1076 +    real (kind=dp), intent(in) :: r, mu
1077 +    real (kind=dp), intent(inout) :: scale, dscale
1078 +    real (kind=dp) :: rl, ru, mulow, minRatio, temp, scaleVal
1079 +
1080 +    ! distances must be in angstroms
1081 +    rl = 2.75d0
1082 +    ru = 3.75d0
1083 +    mulow = 0.0d0 !3.3856d0 ! 1.84 * 1.84
1084 +    minRatio = mulow / (mu*mu)
1085 +    scaleVal = 1.0d0 - minRatio
1086 +    
1087 +    if (r.lt.rl) then
1088 +       scale = minRatio
1089 +       dscale = 0.0d0
1090 +    elseif (r.gt.ru) then
1091 +       scale = 1.0d0
1092 +       dscale = 0.0d0
1093 +    else
1094 +       scale = 1.0d0 - scaleVal*((ru + 2.0d0*r - 3.0d0*rl) * (ru-r)**2) &
1095 +                        / ((ru - rl)**3)
1096 +       dscale = -scaleVal * 6.0d0 * (r-ru)*(r-rl)/((ru - rl)**3)    
1097 +    endif
1098 +        
1099 +    return
1100 +  end subroutine calc_switch
1101 +
1102    subroutine destroyElectrostaticTypes()
1103  
1104      if(allocated(ElectrostaticMap)) deallocate(ElectrostaticMap)

Diff Legend

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