553 |
|
real(kind=dp) :: Ri, Ri3, Ri6, Ri7, Ri12, Ri13, R126, R137, prefactor |
554 |
|
real(kind=dp) :: chipoalphap2, chioalpha2, ec, epsnot |
555 |
|
real(kind=dp) :: drdotudx, drdotudy, drdotudz |
556 |
+ |
real(kind=dp) :: drdotudux, drdotuduy, drdotuduz |
557 |
|
real(kind=dp) :: ljeps, ljsigma, sigmaratio, sigmaratioi |
558 |
|
integer :: ljt1, ljt2, atid1, atid2, gbt1, gbt2 |
559 |
|
logical :: gb_first |
628 |
|
ljsigma = getSigma(atid1) |
629 |
|
ljeps = getEpsilon(atid1) |
630 |
|
endif |
631 |
< |
|
631 |
> |
|
632 |
> |
write(*,*) 'd u', dx, dy, dz, ul(1), ul(2), ul(3) |
633 |
> |
|
634 |
|
rdotu = (dx*ul(1)+dy*ul(2)+dz*ul(3))*ri |
635 |
< |
|
635 |
> |
|
636 |
|
drdotudx = ul(1)*ri-rdotu*dx*ri*ri |
637 |
|
drdotudy = ul(2)*ri-rdotu*dy*ri*ri |
638 |
|
drdotudz = ul(3)*ri-rdotu*dz*ri*ri |
639 |
+ |
drdotudux = drdx |
640 |
+ |
drdotuduy = drdy |
641 |
+ |
drdotuduz = drdz |
642 |
+ |
|
643 |
|
|
644 |
|
moom = 1.0d0 / gb_mu |
645 |
|
mum1 = gb_mu-1 |
669 |
|
|
670 |
|
mess = 1-rdotu*rdotu*chioalpha2 |
671 |
|
sab = 1.0d0/dsqrt(mess) |
672 |
+ |
|
673 |
+ |
write(*,*) 's', sc, sab, rdotu, chioalpha2 |
674 |
|
dsabdct = sc*sab*sab*sab*rdotu*chioalpha2 |
675 |
|
|
676 |
|
eab = 1-chipoalphap2*rdotu*rdotu |
677 |
|
eabf = enot*eab**gb_mu |
678 |
+ |
|
679 |
+ |
write(*,*) 'e', enot, chipoalphap2, gb_mu, rdotu, eab, mum1 |
680 |
+ |
|
681 |
|
depmudct = -2*enot*chipoalphap2*gb_mu*rdotu*eab**mum1 |
682 |
|
|
683 |
|
BigR = (r - sab*sc + sc)/sc |
684 |
|
dBigRdx = (drdx -dsabdct*drdotudx)/sc |
685 |
|
dBigRdy = (drdy -dsabdct*drdotudy)/sc |
686 |
|
dBigRdz = (drdz -dsabdct*drdotudz)/sc |
687 |
< |
dBigRdux = (-dsabdct*drdx)/sc |
688 |
< |
dBigRduy = (-dsabdct*drdy)/sc |
689 |
< |
dBigRduz = (-dsabdct*drdz)/sc |
687 |
> |
dBigRdux = (-dsabdct*drdotudux)/sc |
688 |
> |
dBigRduy = (-dsabdct*drdotuduy)/sc |
689 |
> |
dBigRduz = (-dsabdct*drdotuduz)/sc |
690 |
|
|
691 |
+ |
write(*,*) 'ds dep', dsabdct, depmudct |
692 |
+ |
write(*,*) 'drdotudu', drdotudux, drdotuduy, drdotuduz |
693 |
|
depmudx = depmudct*drdotudx |
694 |
|
depmudy = depmudct*drdotudy |
695 |
|
depmudz = depmudct*drdotudz |
696 |
< |
depmudux = depmudct*drdx |
697 |
< |
depmuduy = depmudct*drdy |
698 |
< |
depmuduz = depmudct*drdz |
696 |
> |
depmudux = depmudct*drdotudux |
697 |
> |
depmuduy = depmudct*drdotuduy |
698 |
> |
depmuduz = depmudct*drdotuduz |
699 |
|
|
700 |
|
Ri = 1.0d0/BigR |
701 |
|
Ri3 = Ri*Ri*Ri |
711 |
|
dUdx = prefactor*(eabf*R137*dBigRdx + R126*depmudx) |
712 |
|
dUdy = prefactor*(eabf*R137*dBigRdy + R126*depmudy) |
713 |
|
dUdz = prefactor*(eabf*R137*dBigRdz + R126*depmudz) |
714 |
< |
write(*,*) 'p', prefactor, eabf, r137,dbigrdux, depmudux, r126 |
714 |
> |
write(*,*) 'dRdu', dbigrdux, dbigrduy, dbigrduz |
715 |
> |
write(*,*) 'dEdu', depmudux, depmuduy, depmuduz |
716 |
|
dUdux = prefactor*(eabf*R137*dBigRdux + R126*depmudux) |
717 |
|
dUduy = prefactor*(eabf*R137*dBigRduy + R126*depmuduy) |
718 |
|
dUduz = prefactor*(eabf*R137*dBigRduz + R126*depmuduz) |
752 |
|
t(1,atom1) = t(1,atom1) + ul(2)*dUduz - ul(3)*dUduy |
753 |
|
t(2,atom1) = t(2,atom1) + ul(3)*dUdux - ul(1)*dUduz |
754 |
|
t(3,atom1) = t(3,atom1) + ul(1)*dUduy - ul(2)*dUdux |
755 |
< |
write(*,*) t(1,atom1), t(2,atom1), t(3,atom1) |
755 |
> |
write(*,*) 'T1', t(1,atom1), t(2,atom1), t(3,atom1) |
756 |
|
else |
757 |
|
t(1,atom2) = t(1,atom2) + ul(2)*dUduz - ul(3)*dUduy |
758 |
|
t(2,atom2) = t(2,atom2) + ul(3)*dUdux - ul(1)*dUduz |
759 |
|
t(3,atom2) = t(3,atom2) + ul(1)*dUduy - ul(2)*dUdux |
760 |
|
|
761 |
< |
write(*,*) t(1,atom2), t(2,atom2), t(3,atom2) |
761 |
> |
write(*,*) 'T2', t(1,atom2), t(2,atom2), t(3,atom2) |
762 |
|
endif |
763 |
|
|
764 |
|
#endif |