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 2375 by gezelter, Mon Oct 17 19:12:45 2005 UTC vs.
Revision 2379 by kdaily, Mon Oct 17 22:42:07 2005 UTC

# Line 539 | Line 539 | contains
539      logical, intent(in) :: do_pot
540      real (kind = dp), dimension(3) :: ul
541      
542 <    real(kind=dp) :: gb_sigma, gb_eps, gb_eps_ratio, gb_mu, gb_sigma_l
543 <    real(kind=dp) :: spar, sperp, slj, par2, perp2, sc, slj2
544 <    real(kind=dp) :: s_par2mperp2, s_lj2ppar2
545 <    real(kind=dp) :: enot, eperp, epar, eab, eabf,moom, mum1
546 <    real(kind=dp) :: dx, dy, dz, drdx,drdy,drdz, rdotu
547 <    real(kind=dp) :: mess, sab, dsabdct, eprime, deprimedct, depmudct
542 >    real(kind=dp) :: gb_sigma, gb_eps, gb_eps_ratio, gb_mu, gb_l2b_ratio
543 >    real(kind=dp) :: s0, l2, d2, lj2
544 >    real(kind=dp) :: eE, eS, eab, eabf, moom, mum1
545 >    real(kind=dp) :: dx, dy, dz, drdx, drdy, drdz, rdotu
546 >    real(kind=dp) :: mess, sab, dsabdct, depmudct
547      real(kind=dp) :: epmu, depmudx, depmudy, depmudz
548      real(kind=dp) :: depmudux, depmuduy, depmuduz
549      real(kind=dp) :: BigR, dBigRdx, dBigRdy, dBigRdz
# Line 553 | Line 552 | contains
552      real(kind=dp) :: Ri, Ri3, Ri6, Ri7, Ri12, Ri13, R126, R137, prefactor
553      real(kind=dp) :: chipoalphap2, chioalpha2, ec, epsnot
554      real(kind=dp) :: drdotudx, drdotudy, drdotudz    
555 <    real(kind=dp) :: ljeps, ljsigma, sigmaratio, sigmaratioi
555 >    real(kind=dp) :: drdotudux, drdotuduy, drdotuduz    
556 >    real(kind=dp) :: ljeps, ljsigma
557      integer :: ljt1, ljt2, atid1, atid2, gbt1, gbt2
558      logical :: gb_first
559      
# Line 601 | Line 601 | contains
601         ul(3) = A(9,atom1)      
602   #endif
603         gb_sigma     = GBMap%GBtypes(gbt1)%sigma      
604 +       gb_l2b_ratio = GBMap%GBtypes(gbt1)%l2b_ratio
605         gb_eps       = GBMap%GBtypes(gbt1)%eps        
606         gb_eps_ratio = GBMap%GBtypes(gbt1)%eps_ratio  
607         gb_mu        = GBMap%GBtypes(gbt1)%mu        
607 –       gb_sigma_l   = GBMap%GBtypes(gbt1)%sigma_l
608  
609         ljsigma = getSigma(atid2)
610         ljeps = getEpsilon(atid2)
# Line 619 | Line 619 | contains
619         ul(3) = A(9,atom2)      
620   #endif
621         gb_sigma     = GBMap%GBtypes(gbt2)%sigma      
622 +       gb_l2b_ratio = GBMap%GBtypes(gbt2)%l2b_ratio
623         gb_eps       = GBMap%GBtypes(gbt2)%eps        
624         gb_eps_ratio = GBMap%GBtypes(gbt2)%eps_ratio  
625         gb_mu        = GBMap%GBtypes(gbt2)%mu        
625 –       gb_sigma_l   = GBMap%GBtypes(gbt2)%sigma_l
626  
627         ljsigma = getSigma(atid1)
628         ljeps = getEpsilon(atid1)
629      endif
630 <    
630 >  
631 >    write(*,*) 'd u', dx, dy, dz, ul(1), ul(2), ul(3)
632 >
633      rdotu = (dx*ul(1)+dy*ul(2)+dz*ul(3))*ri
634 <    
634 >  
635      drdotudx = ul(1)*ri-rdotu*dx*ri*ri
636      drdotudy = ul(2)*ri-rdotu*dy*ri*ri
637      drdotudz = ul(3)*ri-rdotu*dz*ri*ri
638 <    
639 <    moom =  1.0d0 / gb_mu
640 <    mum1 = gb_mu-1
639 <    
640 <    sperp = gb_sigma
641 <    spar =  gb_sigma_l
642 <    slj = ljsigma
643 <    slj2 = slj*slj
638 >    drdotudux = drdx
639 >    drdotuduy = drdy
640 >    drdotuduz = drdz
641  
642 <    chioalpha2 =1-((sperp+slj)*(sperp+slj))/((spar+slj)*(spar+slj))
643 <    sc = (sperp+slj)/2.0d0
644 <  
645 <    par2 = spar*spar
649 <    perp2 = sperp*sperp
650 <    s_par2mperp2 = par2 - perp2
651 <    s_lj2ppar2 = slj2 + par2
642 >    l2 = (gb_sigma*gb_l2b_ratio)**2
643 >    d2 = gb_sigma**2
644 >    lj2 = ljsigma**2
645 >    s0 = dsqrt(d2 + lj2)
646  
647 <    !! check these ! left from old code
654 <    !! kdaily e0 may need to be (gb_eps + lj_eps)/2
655 <    
656 <    eperp = dsqrt(gb_eps*ljeps)
657 <    epar = eperp*gb_eps_ratio
658 <    enot = dsqrt(ljeps*eperp)
659 <    chipoalphap2 = 1+(dsqrt(epar*ljeps)/enot)**moom
647 >    chioalpha2 = (l2 - d2)/(l2 + lj2)
648  
649 +    eE = dsqrt(gb_eps*ljeps)
650 +    eS = dsqrt(gb_eps*gb_eps_ratio*ljeps)
651 +    moom =  1.0d0 / gb_mu
652 +    mum1 = gb_mu-1
653 +    chipoalphap2 = 1 - (eE/eS)**moom
654 +
655      !! mess matches cleaver (eq 20)
656      
657      mess = 1-rdotu*rdotu*chioalpha2
658      sab = 1.0d0/dsqrt(mess)
659 <    dsabdct = sc*sab*sab*sab*rdotu*chioalpha2
659 >
660 >    write(*,*) 's', s0, sab, rdotu, chioalpha2
661 >    dsabdct = s0*sab*sab*sab*rdotu*chioalpha2
662        
663      eab = 1-chipoalphap2*rdotu*rdotu
664 <    eabf = enot*eab**gb_mu
665 <    depmudct = -2*enot*chipoalphap2*gb_mu*rdotu*eab**mum1
664 >    eabf = eS*(eab**gb_mu)
665 >
666 >    write(*,*)  'e', eS, chipoalphap2, gb_mu, rdotu, eab, mum1
667 >
668 >    depmudct = -2*eS*chipoalphap2*gb_mu*rdotu*(eab**mum1)
669          
670 <    BigR = (r - sab*sc + sc)/sc
671 <    dBigRdx = (drdx -dsabdct*drdotudx)/sc
672 <    dBigRdy = (drdy -dsabdct*drdotudy)/sc
673 <    dBigRdz = (drdz -dsabdct*drdotudz)/sc
674 <    dBigRdux = (-dsabdct*drdx)/sc
675 <    dBigRduy = (-dsabdct*drdy)/sc
676 <    dBigRduz = (-dsabdct*drdz)/sc
670 >    BigR = (r - sab*s0 + s0)/s0
671 >    dBigRdx = (drdx -dsabdct*drdotudx)/s0
672 >    dBigRdy = (drdy -dsabdct*drdotudy)/s0
673 >    dBigRdz = (drdz -dsabdct*drdotudz)/s0
674 >    dBigRdux = (-dsabdct*drdotudux)/s0
675 >    dBigRduy = (-dsabdct*drdotuduy)/s0
676 >    dBigRduz = (-dsabdct*drdotuduz)/s0
677      
678 +    write(*,*) 'ds dep', dsabdct, depmudct
679 +    write(*,*) 'drdotudu', drdotudux, drdotuduy, drdotuduz
680      depmudx = depmudct*drdotudx
681      depmudy = depmudct*drdotudy
682      depmudz = depmudct*drdotudz
683 <    depmudux = depmudct*drdx
684 <    depmuduy = depmudct*drdy
685 <    depmuduz = depmudct*drdz
683 >    depmudux = depmudct*drdotudux
684 >    depmuduy = depmudct*drdotuduy
685 >    depmuduz = depmudct*drdotuduz
686      
687      Ri = 1.0d0/BigR
688      Ri3 = Ri*Ri*Ri
# Line 697 | Line 698 | contains
698      dUdx = prefactor*(eabf*R137*dBigRdx + R126*depmudx)
699      dUdy = prefactor*(eabf*R137*dBigRdy + R126*depmudy)
700      dUdz = prefactor*(eabf*R137*dBigRdz + R126*depmudz)
701 <    write(*,*) 'p', prefactor, eabf, r137,dbigrdux, depmudux, r126
701 >    write(*,*) 'dRdu',  dbigrdux, dbigrduy, dbigrduz
702 >    write(*,*) 'dEdu',  depmudux, depmuduy, depmuduz
703      dUdux = prefactor*(eabf*R137*dBigRdux + R126*depmudux)
704      dUduy = prefactor*(eabf*R137*dBigRduy + R126*depmuduy)
705      dUduz = prefactor*(eabf*R137*dBigRduz + R126*depmuduz)
# Line 737 | Line 739 | contains
739         t(1,atom1) = t(1,atom1) + ul(2)*dUduz - ul(3)*dUduy
740         t(2,atom1) = t(2,atom1) + ul(3)*dUdux - ul(1)*dUduz
741         t(3,atom1) = t(3,atom1) + ul(1)*dUduy - ul(2)*dUdux
742 <       write(*,*) t(1,atom1), t(2,atom1), t(3,atom1)
742 >       write(*,*) 'T1', t(1,atom1), t(2,atom1), t(3,atom1)
743      else
744         t(1,atom2) = t(1,atom2) + ul(2)*dUduz - ul(3)*dUduy
745         t(2,atom2) = t(2,atom2) + ul(3)*dUdux - ul(1)*dUduz
746         t(3,atom2) = t(3,atom2) + ul(1)*dUduy - ul(2)*dUdux
747  
748 <       write(*,*) t(1,atom2), t(2,atom2), t(3,atom2)
748 >       write(*,*) 'T2', t(1,atom2), t(2,atom2), t(3,atom2)
749      endif
750  
751   #endif

Diff Legend

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