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 |
|
|
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) |
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 |
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 |
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 |
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 |
|
|
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 |
|
|
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 = dsqrt(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 = 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 |
650 |
|
|
651 |
|
!! mess matches cleaver (eq 20) |
652 |
|
|
653 |
|
mess = 1-rdotu*rdotu*chioalpha2 |
654 |
|
sab = 1.0d0/dsqrt(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 |
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 |
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) |