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

Comparing trunk/OOPSE-3.0/src/UseTheForce/DarkSide/gb.F90 (file contents):
Revision 2379 by kdaily, Mon Oct 17 22:42:07 2005 UTC vs.
Revision 2385 by gezelter, Tue Oct 18 21:13:49 2005 UTC

# 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 418 | Line 418 | contains
418      gmum = gmu*gpi
419  
420      curlyE = 1.0d0/dsqrt(1.0d0 - Chi*Chi*u1dotu2*u1dotu2)
421 – !!$
422 – !!$    dcE = -(curlyE**3)*Chi*Chi*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 592 | 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      
# Line 610 | Line 608 | contains
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
# Line 626 | Line 624 | contains
624  
625         ljsigma = getSigma(atid1)
626         ljeps = getEpsilon(atid1)
627 <    endif
630 <  
631 <    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 646 | Line 642 | contains
642  
643      chioalpha2 = (l2 - d2)/(l2 + lj2)
644  
645 <    eE = dsqrt(gb_eps*ljeps)
646 <    eS = dsqrt(gb_eps*gb_eps_ratio*ljeps)
645 >    eE = dsqrt(gb_eps*gb_eps_ratio*ljeps)
646 >    eS = dsqrt(gb_eps*ljeps)
647      moom =  1.0d0 / gb_mu
648      mum1 = gb_mu-1
649      chipoalphap2 = 1 - (eE/eS)**moom
# Line 657 | Line 653 | contains
653      mess = 1-rdotu*rdotu*chioalpha2
654      sab = 1.0d0/dsqrt(mess)
655  
660 –    write(*,*) 's', s0, sab, rdotu, chioalpha2
656      dsabdct = s0*sab*sab*sab*rdotu*chioalpha2
657        
658      eab = 1-chipoalphap2*rdotu*rdotu
659      eabf = eS*(eab**gb_mu)
660  
666 –    write(*,*)  'e', eS, chipoalphap2, gb_mu, rdotu, eab, mum1
667 –
661      depmudct = -2*eS*chipoalphap2*gb_mu*rdotu*(eab**mum1)
662          
663      BigR = (r - sab*s0 + s0)/s0
# Line 675 | Line 668 | contains
668      dBigRduy = (-dsabdct*drdotuduy)/s0
669      dBigRduz = (-dsabdct*drdotuduz)/s0
670      
678 –    write(*,*) 'ds dep', dsabdct, depmudct
679 –    write(*,*) 'drdotudu', drdotudux, drdotuduy, drdotuduz
671      depmudx = depmudct*drdotudx
672      depmudy = depmudct*drdotudy
673      depmudz = depmudct*drdotudz
# Line 695 | 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)
705 <    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 731 | 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:
735 <    
736 <    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
742 <       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
747 <
748 <       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

Diff Legend

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