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 |
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 |
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 |
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 |
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 |
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 |
|
|
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 |
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 |
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 |
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 |
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 |
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) |