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

Comparing trunk/OOPSE-4/src/UseTheForce/DarkSide/gb.F90 (file contents):
Revision 2378 by gezelter, Mon Oct 17 21:47:45 2005 UTC vs.
Revision 2756 by gezelter, Wed May 17 15:37:15 2006 UTC

# Line 229 | Line 229 | contains
229         sigma_l   = GBMap%GBtypes(gbt1)%sigma_l    
230         eps_l     = GBMap%GBtypes(gbt1)%eps_l      
231      else
232 <       call handleError("GB", "GB-pair was called with two different GB types! OOPSE can only handle interactions for one GB type at a time.")
232 >       call handleError("GB", "GB-pair was called with two different GB types!")
233      endif
234  
235      s2 = (l2b_ratio)**2
# Line 241 | Line 241 | contains
241      r4 = r2*r2
242  
243   #ifdef IS_MPI
244 <    ul1(1) = A_Row(3,atom1)
245 <    ul1(2) = A_Row(6,atom1)
244 >    ul1(1) = A_Row(7,atom1)
245 >    ul1(2) = A_Row(8,atom1)
246      ul1(3) = A_Row(9,atom1)
247  
248 <    ul2(1) = A_Col(3,atom2)
249 <    ul2(2) = A_Col(6,atom2)
248 >    ul2(1) = A_Col(7,atom2)
249 >    ul2(2) = A_Col(8,atom2)
250      ul2(3) = A_Col(9,atom2)
251   #else
252 <    ul1(1) = A(3,atom1)
253 <    ul1(2) = A(6,atom1)
252 >    ul1(1) = A(7,atom1)
253 >    ul1(2) = A(8,atom1)
254      ul1(3) = A(9,atom1)
255  
256 <    ul2(1) = A(3,atom2)
257 <    ul2(2) = A(6,atom2)
256 >    ul2(1) = A(7,atom2)
257 >    ul2(2) = A(8,atom2)
258      ul2(3) = A(9,atom2)
259   #endif
260      
# Line 417 | Line 417 | contains
417      gpi = 1.0d0 / gp
418      gmum = gmu*gpi
419  
420 <    curlyE = 1.0d0/dsqrt(1.0d0 - Chi*Chi*u1dotu2*u1dotu2)
421 < !!$
422 < !!$    dcE = -(curlyE**3)*Chi*Chi*u1dotu2
420 >    curlyE = 1.0d0/sqrt(1.0d0 - Chi*Chi*u1dotu2*u1dotu2)
421      dcE = (curlyE**3)*Chi*Chi*u1dotu2
422  
423      dcEdu1x = dcE*ul2(1)
# Line 470 | Line 468 | contains
468      f_Col(2,atom2) = f_Col(2,atom2) - dUdy
469      f_Col(3,atom2) = f_Col(3,atom2) - dUdz
470      
471 <    t_Row(1,atom1) = t_Row(1,atom1)- ul1(3)*dUdu1y + ul1(2)*dUdu1z
472 <    t_Row(2,atom1) = t_Row(2,atom1)- ul1(1)*dUdu1z + ul1(3)*dUdu1x
473 <    t_Row(3,atom1) = t_Row(3,atom1)- ul1(2)*dUdu1x + ul1(1)*dUdu1y
471 >    t_Row(1,atom1) = t_Row(1,atom1) + ul1(3)*dUdu1y - ul1(2)*dUdu1z
472 >    t_Row(2,atom1) = t_Row(2,atom1) + ul1(1)*dUdu1z - ul1(3)*dUdu1x
473 >    t_Row(3,atom1) = t_Row(3,atom1) + ul1(2)*dUdu1x - ul1(1)*dUdu1y
474      
475 <    t_Col(1,atom2) = t_Col(1,atom2) - ul2(3)*dUdu2y + ul2(2)*dUdu2z
476 <    t_Col(2,atom2) = t_Col(2,atom2) - ul2(1)*dUdu2z + ul2(3)*dUdu2x
477 <    t_Col(3,atom2) = t_Col(3,atom2) - ul2(2)*dUdu2x + ul2(1)*dUdu2y
475 >    t_Col(1,atom2) = t_Col(1,atom2) + ul2(3)*dUdu2y - ul2(2)*dUdu2z
476 >    t_Col(2,atom2) = t_Col(2,atom2) + ul2(1)*dUdu2z - ul2(3)*dUdu2x
477 >    t_Col(3,atom2) = t_Col(3,atom2) + ul2(2)*dUdu2x - ul2(1)*dUdu2y
478   #else
479      f(1,atom1) = f(1,atom1) + dUdx
480      f(2,atom1) = f(2,atom1) + dUdy
# Line 486 | Line 484 | contains
484      f(2,atom2) = f(2,atom2) - dUdy
485      f(3,atom2) = f(3,atom2) - dUdz
486      
487 <    t(1,atom1) = t(1,atom1)- ul1(3)*dUdu1y + ul1(2)*dUdu1z
488 <    t(2,atom1) = t(2,atom1)- ul1(1)*dUdu1z + ul1(3)*dUdu1x
489 <    t(3,atom1) = t(3,atom1)- ul1(2)*dUdu1x + ul1(1)*dUdu1y
487 >    t(1,atom1) = t(1,atom1) + ul1(3)*dUdu1y - ul1(2)*dUdu1z
488 >    t(2,atom1) = t(2,atom1) + ul1(1)*dUdu1z - ul1(3)*dUdu1x
489 >    t(3,atom1) = t(3,atom1) + ul1(2)*dUdu1x - ul1(1)*dUdu1y
490      
491 <    t(1,atom2) = t(1,atom2)- ul2(3)*dUdu2y + ul2(2)*dUdu2z
492 <    t(2,atom2) = t(2,atom2)- ul2(1)*dUdu2z + ul2(3)*dUdu2x
493 <    t(3,atom2) = t(3,atom2)- ul2(2)*dUdu2x + ul2(1)*dUdu2y
491 >    t(1,atom2) = t(1,atom2) + ul2(3)*dUdu2y - ul2(2)*dUdu2z
492 >    t(2,atom2) = t(2,atom2) + ul2(1)*dUdu2z - ul2(3)*dUdu2x
493 >    t(3,atom2) = t(3,atom2) + ul2(2)*dUdu2x - ul2(1)*dUdu2y
494   #endif
495    
496      if (do_pot) then
# Line 524 | Line 522 | contains
522      return
523    end subroutine do_gb_pair
524  
525 <  subroutine do_gb_lj_pair(atom1, atom2, d, r, r2, sw, vpair, fpair, &
525 >  subroutine do_gb_lj_pair(atom1, atom2, d, r, r2, rcut, sw, vpair, fpair, &
526         pot, A, f, t, do_pot)
527      
528      integer, intent(in) :: atom1, atom2
529      integer :: id1, id2
530 <    real (kind=dp), intent(inout) :: r, r2
530 >    real (kind=dp), intent(inout) :: r, r2, rcut
531      real (kind=dp), dimension(3), intent(in) :: d
532      real (kind=dp), dimension(3), intent(inout) :: fpair
533      real (kind=dp) :: pot, sw, vpair
# Line 539 | Line 537 | contains
537      logical, intent(in) :: do_pot
538      real (kind = dp), dimension(3) :: ul
539      
540 <    real(kind=dp) :: gb_sigma, gb_eps, gb_eps_ratio, gb_mu, gb_sigma_l
541 <    real(kind=dp) :: spar, sperp, slj, par2, perp2, sc, slj2
542 <    real(kind=dp) :: s_par2mperp2, s_lj2ppar2
543 <    real(kind=dp) :: enot, eperp, epar, eab, eabf,moom, mum1
544 <    real(kind=dp) :: dx, dy, dz, drdx,drdy,drdz, rdotu
547 <    real(kind=dp) :: mess, sab, dsabdct, eprime, deprimedct, depmudct
540 >    real(kind=dp) :: gb_sigma, gb_eps, gb_eps_ratio, gb_mu, gb_l2b_ratio
541 >    real(kind=dp) :: s0, l2, d2, lj2
542 >    real(kind=dp) :: eE, eS, eab, eabf, moom, mum1
543 >    real(kind=dp) :: dx, dy, dz, drdx, drdy, drdz, rdotu
544 >    real(kind=dp) :: mess, sab, dsabdct, depmudct
545      real(kind=dp) :: epmu, depmudx, depmudy, depmudz
546      real(kind=dp) :: depmudux, depmuduy, depmuduz
547      real(kind=dp) :: BigR, dBigRdx, dBigRdy, dBigRdz
# Line 554 | Line 551 | contains
551      real(kind=dp) :: chipoalphap2, chioalpha2, ec, epsnot
552      real(kind=dp) :: drdotudx, drdotudy, drdotudz    
553      real(kind=dp) :: drdotudux, drdotuduy, drdotuduz    
554 <    real(kind=dp) :: ljeps, ljsigma, sigmaratio, sigmaratioi
554 >    real(kind=dp) :: ljeps, ljsigma
555      integer :: ljt1, ljt2, atid1, atid2, gbt1, gbt2
556      logical :: gb_first
557      
# Line 593 | Line 590 | contains
590      
591      if(gb_first)then
592   #ifdef IS_MPI
593 <       ul(1) = A_Row(3,atom1)
594 <       ul(2) = A_Row(6,atom1)
593 >       ul(1) = A_Row(7,atom1)
594 >       ul(2) = A_Row(8,atom1)
595         ul(3) = A_Row(9,atom1)
596   #else
597 <       ul(1) = A(3,atom1)
598 <       ul(2) = A(6,atom1)
597 >       ul(1) = A(7,atom1)
598 >       ul(2) = A(8,atom1)
599         ul(3) = A(9,atom1)      
600   #endif
601         gb_sigma     = GBMap%GBtypes(gbt1)%sigma      
602 +       gb_l2b_ratio = GBMap%GBtypes(gbt1)%l2b_ratio
603         gb_eps       = GBMap%GBtypes(gbt1)%eps        
604         gb_eps_ratio = GBMap%GBtypes(gbt1)%eps_ratio  
605         gb_mu        = GBMap%GBtypes(gbt1)%mu        
608 –       gb_sigma_l   = GBMap%GBtypes(gbt1)%sigma_l
606  
607         ljsigma = getSigma(atid2)
608         ljeps = getEpsilon(atid2)
609      else
610   #ifdef IS_MPI
611 <       ul(1) = A_Col(3,atom2)
612 <       ul(2) = A_Col(6,atom2)
611 >       ul(1) = A_Col(7,atom2)
612 >       ul(2) = A_Col(8,atom2)
613         ul(3) = A_Col(9,atom2)
614   #else
615 <       ul(1) = A(3,atom2)
616 <       ul(2) = A(6,atom2)
617 <       ul(3) = A(9,atom2)      
615 >       ul(1) = A(7,atom2)
616 >       ul(2) = A(8,atom2)
617 >       ul(3) = A(9,atom2)    
618   #endif
619         gb_sigma     = GBMap%GBtypes(gbt2)%sigma      
620 +       gb_l2b_ratio = GBMap%GBtypes(gbt2)%l2b_ratio
621         gb_eps       = GBMap%GBtypes(gbt2)%eps        
622         gb_eps_ratio = GBMap%GBtypes(gbt2)%eps_ratio  
623         gb_mu        = GBMap%GBtypes(gbt2)%mu        
626 –       gb_sigma_l   = GBMap%GBtypes(gbt2)%sigma_l
624  
625         ljsigma = getSigma(atid1)
626         ljeps = getEpsilon(atid1)
627 <    endif
631 <  
632 <    write(*,*) 'd u', dx, dy, dz, ul(1), ul(2), ul(3)
627 >    endif  
628  
629      rdotu = (dx*ul(1)+dy*ul(2)+dz*ul(3))*ri
630    
# Line 640 | Line 635 | contains
635      drdotuduy = drdy
636      drdotuduz = drdz
637  
638 <    
639 <    moom =  1.0d0 / gb_mu
640 <    mum1 = gb_mu-1
641 <    
647 <    sperp = gb_sigma
648 <    spar =  gb_sigma_l
649 <    slj = ljsigma
650 <    slj2 = slj*slj
638 >    l2 = (gb_sigma*gb_l2b_ratio)**2
639 >    d2 = gb_sigma**2
640 >    lj2 = ljsigma**2
641 >    s0 = sqrt(d2 + lj2)
642  
643 <    chioalpha2 =1-((sperp+slj)*(sperp+slj))/((spar+slj)*(spar+slj))
653 <    sc = (sperp+slj)/2.0d0
654 <  
655 <    par2 = spar*spar
656 <    perp2 = sperp*sperp
657 <    s_par2mperp2 = par2 - perp2
658 <    s_lj2ppar2 = slj2 + par2
643 >    chioalpha2 = (l2 - d2)/(l2 + lj2)
644  
645 <    !! check these ! left from old code
646 <    !! kdaily e0 may need to be (gb_eps + lj_eps)/2
647 <    
648 <    eperp = dsqrt(gb_eps*ljeps)
649 <    epar = eperp*gb_eps_ratio
665 <    enot = dsqrt(ljeps*eperp)
666 <    chipoalphap2 = 1+(dsqrt(epar*ljeps)/enot)**moom
645 >    eE = sqrt(gb_eps*gb_eps_ratio*ljeps)
646 >    eS = sqrt(gb_eps*ljeps)
647 >    moom =  1.0d0 / gb_mu
648 >    mum1 = gb_mu-1
649 >    chipoalphap2 = 1 - (eE/eS)**moom
650  
651      !! mess matches cleaver (eq 20)
652      
653      mess = 1-rdotu*rdotu*chioalpha2
654 <    sab = 1.0d0/dsqrt(mess)
654 >    sab = 1.0d0/sqrt(mess)
655  
656 <    write(*,*) 's', sc, sab, rdotu, chioalpha2
674 <    dsabdct = sc*sab*sab*sab*rdotu*chioalpha2
656 >    dsabdct = s0*sab*sab*sab*rdotu*chioalpha2
657        
658      eab = 1-chipoalphap2*rdotu*rdotu
659 <    eabf = enot*eab**gb_mu
659 >    eabf = eS*(eab**gb_mu)
660  
661 <    write(*,*)  'e', enot, chipoalphap2, gb_mu, rdotu, eab, mum1
680 <
681 <    depmudct = -2*enot*chipoalphap2*gb_mu*rdotu*eab**mum1
661 >    depmudct = -2*eS*chipoalphap2*gb_mu*rdotu*(eab**mum1)
662          
663 <    BigR = (r - sab*sc + sc)/sc
664 <    dBigRdx = (drdx -dsabdct*drdotudx)/sc
665 <    dBigRdy = (drdy -dsabdct*drdotudy)/sc
666 <    dBigRdz = (drdz -dsabdct*drdotudz)/sc
667 <    dBigRdux = (-dsabdct*drdotudux)/sc
668 <    dBigRduy = (-dsabdct*drdotuduy)/sc
669 <    dBigRduz = (-dsabdct*drdotuduz)/sc
663 >    BigR = (r - sab*s0 + s0)/s0
664 >    dBigRdx = (drdx -dsabdct*drdotudx)/s0
665 >    dBigRdy = (drdy -dsabdct*drdotudy)/s0
666 >    dBigRdz = (drdz -dsabdct*drdotudz)/s0
667 >    dBigRdux = (-dsabdct*drdotudux)/s0
668 >    dBigRduy = (-dsabdct*drdotuduy)/s0
669 >    dBigRduz = (-dsabdct*drdotuduz)/s0
670      
691 –    write(*,*) 'ds dep', dsabdct, depmudct
692 –    write(*,*) 'drdotudu', drdotudux, drdotuduy, drdotuduz
671      depmudx = depmudct*drdotudx
672      depmudy = depmudct*drdotudy
673      depmudz = depmudct*drdotudz
# Line 708 | Line 686 | contains
686      
687      prefactor = 4.0d0
688      
689 <    dUdx = prefactor*(eabf*R137*dBigRdx + R126*depmudx)
690 <    dUdy = prefactor*(eabf*R137*dBigRdy + R126*depmudy)
691 <    dUdz = prefactor*(eabf*R137*dBigRdz + R126*depmudz)
692 <    write(*,*) 'dRdu',  dbigrdux, dbigrduy, dbigrduz
693 <    write(*,*) 'dEdu',  depmudux, depmuduy, depmuduz
694 <    dUdux = prefactor*(eabf*R137*dBigRdux + R126*depmudux)
695 <    dUduy = prefactor*(eabf*R137*dBigRduy + R126*depmuduy)
718 <    dUduz = prefactor*(eabf*R137*dBigRduz + R126*depmuduz)
689 >    dUdx = prefactor*(eabf*R137*dBigRdx + R126*depmudx)*sw
690 >    dUdy = prefactor*(eabf*R137*dBigRdy + R126*depmudy)*sw
691 >    dUdz = prefactor*(eabf*R137*dBigRdz + R126*depmudz)*sw
692 >
693 >    dUdux = prefactor*(eabf*R137*dBigRdux + R126*depmudux)*sw
694 >    dUduy = prefactor*(eabf*R137*dBigRduy + R126*depmuduy)*sw
695 >    dUduz = prefactor*(eabf*R137*dBigRduz + R126*depmuduz)*sw
696      
697   #ifdef IS_MPI
698 <    f_Row(1,atom1) = f_Row(1,atom1) - dUdx
699 <    f_Row(2,atom1) = f_Row(2,atom1) - dUdy
700 <    f_Row(3,atom1) = f_Row(3,atom1) - dUdz
698 >    f_Row(1,atom1) = f_Row(1,atom1) + dUdx
699 >    f_Row(2,atom1) = f_Row(2,atom1) + dUdy
700 >    f_Row(3,atom1) = f_Row(3,atom1) + dUdz
701      
702 <    f_Col(1,atom2) = f_Col(1,atom2) + dUdx
703 <    f_Col(2,atom2) = f_Col(2,atom2) + dUdy
704 <    f_Col(3,atom2) = f_Col(3,atom2) + dUdz
702 >    f_Col(1,atom2) = f_Col(1,atom2) - dUdx
703 >    f_Col(2,atom2) = f_Col(2,atom2) - dUdy
704 >    f_Col(3,atom2) = f_Col(3,atom2) - dUdz
705      
706      if (gb_first) then
707 <       t_Row(1,atom1) = t_Row(1,atom1) + ul(2)*dUduz - ul(3)*dUduy
708 <       t_Row(2,atom1) = t_Row(2,atom1) + ul(3)*dUdux - ul(1)*dUduz
709 <       t_Row(3,atom1) = t_Row(3,atom1) + ul(1)*dUduy - ul(2)*dUdux
707 >       t_Row(1,atom1) = t_Row(1,atom1) - ul(2)*dUduz + ul(3)*dUduy
708 >       t_Row(2,atom1) = t_Row(2,atom1) - ul(3)*dUdux + ul(1)*dUduz
709 >       t_Row(3,atom1) = t_Row(3,atom1) - ul(1)*dUduy + ul(2)*dUdux
710      else
711 <       t_Col(1,atom2) = t_Col(1,atom2) + ul(2)*dUduz - ul(3)*dUduy
712 <       t_Col(2,atom2) = t_Col(2,atom2) + ul(3)*dUdux - ul(1)*dUduz
713 <       t_Col(3,atom2) = t_Col(3,atom2) + ul(1)*dUduy - ul(2)*dUdux
711 >       t_Col(1,atom2) = t_Col(1,atom2) - ul(2)*dUduz + ul(3)*dUduy
712 >       t_Col(2,atom2) = t_Col(2,atom2) - ul(3)*dUdux + ul(1)*dUduz
713 >       t_Col(3,atom2) = t_Col(3,atom2) - ul(1)*dUduy + ul(2)*dUdux
714      endif
715   #else    
716      f(1,atom1) = f(1,atom1) + dUdx
# Line 744 | Line 721 | contains
721      f(2,atom2) = f(2,atom2) - dUdy
722      f(3,atom2) = f(3,atom2) - dUdz
723      
724 <    ! torques are cross products:
748 <    
749 <    write(*,*) 'dU', dUdux, duduy, duduz
724 >    ! torques are cross products:    
725  
726      if (gb_first) then
727 <       t(1,atom1) = t(1,atom1) + ul(2)*dUduz - ul(3)*dUduy
728 <       t(2,atom1) = t(2,atom1) + ul(3)*dUdux - ul(1)*dUduz
729 <       t(3,atom1) = t(3,atom1) + ul(1)*dUduy - ul(2)*dUdux
755 <       write(*,*) 'T1', t(1,atom1), t(2,atom1), t(3,atom1)
727 >       t(1,atom1) = t(1,atom1) - ul(2)*dUduz + ul(3)*dUduy
728 >       t(2,atom1) = t(2,atom1) - ul(3)*dUdux + ul(1)*dUduz
729 >       t(3,atom1) = t(3,atom1) - ul(1)*dUduy + ul(2)*dUdux
730      else
731 <       t(1,atom2) = t(1,atom2) + ul(2)*dUduz - ul(3)*dUduy
732 <       t(2,atom2) = t(2,atom2) + ul(3)*dUdux - ul(1)*dUduz
733 <       t(3,atom2) = t(3,atom2) + ul(1)*dUduy - ul(2)*dUdux
760 <
761 <       write(*,*) 'T2', t(1,atom2), t(2,atom2), t(3,atom2)
731 >       t(1,atom2) = t(1,atom2) - ul(2)*dUduz + ul(3)*dUduy
732 >       t(2,atom2) = t(2,atom2) - ul(3)*dUdux + ul(1)*dUduz
733 >       t(3,atom2) = t(3,atom2) - ul(1)*dUduy + ul(2)*dUdux
734      endif
735  
736   #endif
737        
738      if (do_pot) then
739   #ifdef IS_MPI
740 <       pot_row(VDW_POT,atom1) = pot_row(VDW_POT,atom1) + 2.0d0*eps*R126*sw
741 <       pot_col(VDW_POT,atom2) = pot_col(VDW_POT,atom2) + 2.0d0*eps*R126*sw
740 >       pot_row(VDW_POT,atom1) = pot_row(VDW_POT,atom1) + 2.0d0*eabf*R126*sw
741 >       pot_col(VDW_POT,atom2) = pot_col(VDW_POT,atom2) + 2.0d0*eabf*R126*sw
742   #else
743         pot = pot + prefactor*eabf*R126*sw
744   #endif
745      endif
746      
747 <    vpair = vpair + 4.0*epmu*R126
747 >    vpair = vpair + 4.0*eabf*R126
748   #ifdef IS_MPI
749      id1 = AtomRowToGlobal(atom1)
750      id2 = AtomColToGlobal(atom2)

Diff Legend

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