55 |
|
implicit none |
56 |
|
|
57 |
|
PRIVATE |
58 |
< |
|
58 |
> |
#define __FORTRAN90 |
59 |
> |
#include "UseTheForce/DarkSide/fInteractionMap.h" |
60 |
|
INTEGER, PARAMETER:: CHEBYSHEV_TN = 1 |
61 |
|
INTEGER, PARAMETER:: CHEBYSHEV_UN = 2 |
62 |
|
INTEGER, PARAMETER:: LAGUERRE = 3 |
518 |
|
drdxi = -d(1) / rij |
519 |
|
drdyi = -d(2) / rij |
520 |
|
drdzi = -d(3) / rij |
521 |
< |
drduxi = 0.0d0 |
522 |
< |
drduyi = 0.0d0 |
523 |
< |
drduzi = 0.0d0 |
521 |
> |
drduxi = 0.0_dp |
522 |
> |
drduyi = 0.0_dp |
523 |
> |
drduzi = 0.0_dp |
524 |
|
|
525 |
|
drdxj = d(1) / rij |
526 |
|
drdyj = d(2) / rij |
527 |
|
drdzj = d(3) / rij |
528 |
< |
drduxj = 0.0d0 |
529 |
< |
drduyj = 0.0d0 |
530 |
< |
drduzj = 0.0d0 |
528 |
> |
drduxj = 0.0_dp |
529 |
> |
drduyj = 0.0_dp |
530 |
> |
drduzj = 0.0_dp |
531 |
|
|
532 |
|
! find the atom type id (atid) for each atom: |
533 |
|
#ifdef IS_MPI |
549 |
|
sigma_i = ShapeMap%Shapes(st1)%sigma |
550 |
|
s_i = ShapeMap%Shapes(st1)%sigma |
551 |
|
eps_i = ShapeMap%Shapes(st1)%epsilon |
552 |
< |
dsigmaidx = 0.0d0 |
553 |
< |
dsigmaidy = 0.0d0 |
554 |
< |
dsigmaidz = 0.0d0 |
555 |
< |
dsigmaidux = 0.0d0 |
556 |
< |
dsigmaiduy = 0.0d0 |
557 |
< |
dsigmaiduz = 0.0d0 |
558 |
< |
dsidx = 0.0d0 |
559 |
< |
dsidy = 0.0d0 |
560 |
< |
dsidz = 0.0d0 |
561 |
< |
dsidux = 0.0d0 |
562 |
< |
dsiduy = 0.0d0 |
563 |
< |
dsiduz = 0.0d0 |
564 |
< |
depsidx = 0.0d0 |
565 |
< |
depsidy = 0.0d0 |
566 |
< |
depsidz = 0.0d0 |
567 |
< |
depsidux = 0.0d0 |
568 |
< |
depsiduy = 0.0d0 |
569 |
< |
depsiduz = 0.0d0 |
552 |
> |
dsigmaidx = 0.0_dp |
553 |
> |
dsigmaidy = 0.0_dp |
554 |
> |
dsigmaidz = 0.0_dp |
555 |
> |
dsigmaidux = 0.0_dp |
556 |
> |
dsigmaiduy = 0.0_dp |
557 |
> |
dsigmaiduz = 0.0_dp |
558 |
> |
dsidx = 0.0_dp |
559 |
> |
dsidy = 0.0_dp |
560 |
> |
dsidz = 0.0_dp |
561 |
> |
dsidux = 0.0_dp |
562 |
> |
dsiduy = 0.0_dp |
563 |
> |
dsiduz = 0.0_dp |
564 |
> |
depsidx = 0.0_dp |
565 |
> |
depsidy = 0.0_dp |
566 |
> |
depsidz = 0.0_dp |
567 |
> |
depsidux = 0.0_dp |
568 |
> |
depsiduy = 0.0_dp |
569 |
> |
depsiduz = 0.0_dp |
570 |
|
else |
571 |
|
|
572 |
|
#ifdef IS_MPI |
599 |
|
|
600 |
|
dctidx = - zi * xi / r3 |
601 |
|
dctidy = - zi * yi / r3 |
602 |
< |
dctidz = 1.0d0 / rij - zi2 / r3 |
602 |
> |
dctidz = 1.0_dp / rij - zi2 / r3 |
603 |
|
dctidux = yi / rij ! - (zi * xi2) / r3 |
604 |
|
dctiduy = -xi / rij !- (zi * yi2) / r3 |
605 |
< |
dctiduz = 0.0d0 !zi / rij - (zi2 * zi) / r3 |
605 |
> |
dctiduz = 0.0_dp !zi / rij - (zi2 * zi) / r3 |
606 |
|
|
607 |
|
! this is an attempt to try to truncate the singularity when |
608 |
|
! sin(theta) is near 0.0: |
609 |
|
|
610 |
|
sti2 = 1.0_dp - cti*cti |
611 |
< |
if (dabs(sti2) .lt. 1.0d-12) then |
612 |
< |
proji = sqrt(rij * 1.0d-12) |
613 |
< |
dcpidx = 1.0d0 / proji |
614 |
< |
dcpidy = 0.0d0 |
611 |
> |
if (abs(sti2) .lt. 1.0e-12_dp) then |
612 |
> |
proji = sqrt(rij * 1.0e-12_dp) |
613 |
> |
dcpidx = 1.0_dp / proji |
614 |
> |
dcpidy = 0.0_dp |
615 |
|
dcpidux = xi / proji |
616 |
< |
dcpiduy = 0.0d0 |
617 |
< |
dspidx = 0.0d0 |
618 |
< |
dspidy = 1.0d0 / proji |
619 |
< |
dspidux = 0.0d0 |
616 |
> |
dcpiduy = 0.0_dp |
617 |
> |
dspidx = 0.0_dp |
618 |
> |
dspidy = 1.0_dp / proji |
619 |
> |
dspidux = 0.0_dp |
620 |
|
dspiduy = yi / proji |
621 |
|
else |
622 |
|
proji = sqrt(xi2 + yi2) |
623 |
|
proji3 = proji*proji*proji |
624 |
< |
dcpidx = 1.0d0 / proji - xi2 / proji3 |
624 |
> |
dcpidx = 1.0_dp / proji - xi2 / proji3 |
625 |
|
dcpidy = - xi * yi / proji3 |
626 |
|
dcpidux = xi / proji - (xi2 * xi) / proji3 |
627 |
|
dcpiduy = - (xi * yi2) / proji3 |
628 |
|
dspidx = - xi * yi / proji3 |
629 |
< |
dspidy = 1.0d0 / proji - yi2 / proji3 |
629 |
> |
dspidy = 1.0_dp / proji - yi2 / proji3 |
630 |
|
dspidux = - (yi * xi2) / proji3 |
631 |
|
dspiduy = yi / proji - (yi2 * yi) / proji3 |
632 |
|
endif |
633 |
|
|
634 |
|
cpi = xi / proji |
635 |
< |
dcpidz = 0.0d0 |
636 |
< |
dcpiduz = 0.0d0 |
635 |
> |
dcpidz = 0.0_dp |
636 |
> |
dcpiduz = 0.0_dp |
637 |
|
|
638 |
|
spi = yi / proji |
639 |
< |
dspidz = 0.0d0 |
640 |
< |
dspiduz = 0.0d0 |
639 |
> |
dspidz = 0.0_dp |
640 |
> |
dspiduz = 0.0_dp |
641 |
|
|
642 |
|
call Associated_Legendre(cti, ShapeMap%Shapes(st1)%bigM, & |
643 |
|
ShapeMap%Shapes(st1)%bigL, LMAX, & |
648 |
|
call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, MMAX, & |
649 |
|
CHEBYSHEV_UN, um_i, dum_i) |
650 |
|
|
651 |
< |
sigma_i = 0.0d0 |
652 |
< |
s_i = 0.0d0 |
653 |
< |
eps_i = 0.0d0 |
654 |
< |
dsigmaidx = 0.0d0 |
655 |
< |
dsigmaidy = 0.0d0 |
656 |
< |
dsigmaidz = 0.0d0 |
657 |
< |
dsigmaidux = 0.0d0 |
658 |
< |
dsigmaiduy = 0.0d0 |
659 |
< |
dsigmaiduz = 0.0d0 |
660 |
< |
dsidx = 0.0d0 |
661 |
< |
dsidy = 0.0d0 |
662 |
< |
dsidz = 0.0d0 |
663 |
< |
dsidux = 0.0d0 |
664 |
< |
dsiduy = 0.0d0 |
665 |
< |
dsiduz = 0.0d0 |
666 |
< |
depsidx = 0.0d0 |
667 |
< |
depsidy = 0.0d0 |
668 |
< |
depsidz = 0.0d0 |
669 |
< |
depsidux = 0.0d0 |
670 |
< |
depsiduy = 0.0d0 |
671 |
< |
depsiduz = 0.0d0 |
651 |
> |
sigma_i = 0.0_dp |
652 |
> |
s_i = 0.0_dp |
653 |
> |
eps_i = 0.0_dp |
654 |
> |
dsigmaidx = 0.0_dp |
655 |
> |
dsigmaidy = 0.0_dp |
656 |
> |
dsigmaidz = 0.0_dp |
657 |
> |
dsigmaidux = 0.0_dp |
658 |
> |
dsigmaiduy = 0.0_dp |
659 |
> |
dsigmaiduz = 0.0_dp |
660 |
> |
dsidx = 0.0_dp |
661 |
> |
dsidy = 0.0_dp |
662 |
> |
dsidz = 0.0_dp |
663 |
> |
dsidux = 0.0_dp |
664 |
> |
dsiduy = 0.0_dp |
665 |
> |
dsiduz = 0.0_dp |
666 |
> |
depsidx = 0.0_dp |
667 |
> |
depsidy = 0.0_dp |
668 |
> |
depsidz = 0.0_dp |
669 |
> |
depsidux = 0.0_dp |
670 |
> |
depsiduy = 0.0_dp |
671 |
> |
depsiduz = 0.0_dp |
672 |
|
|
673 |
|
do lm = 1, ShapeMap%Shapes(st1)%nContactFuncs |
674 |
|
l = ShapeMap%Shapes(st1)%ContactFuncLValue(lm) |
806 |
|
sigma_j = ShapeMap%Shapes(st2)%sigma |
807 |
|
s_j = ShapeMap%Shapes(st2)%sigma |
808 |
|
eps_j = ShapeMap%Shapes(st2)%epsilon |
809 |
< |
dsigmajdx = 0.0d0 |
810 |
< |
dsigmajdy = 0.0d0 |
811 |
< |
dsigmajdz = 0.0d0 |
812 |
< |
dsigmajdux = 0.0d0 |
813 |
< |
dsigmajduy = 0.0d0 |
814 |
< |
dsigmajduz = 0.0d0 |
815 |
< |
dsjdx = 0.0d0 |
816 |
< |
dsjdy = 0.0d0 |
817 |
< |
dsjdz = 0.0d0 |
818 |
< |
dsjdux = 0.0d0 |
819 |
< |
dsjduy = 0.0d0 |
820 |
< |
dsjduz = 0.0d0 |
821 |
< |
depsjdx = 0.0d0 |
822 |
< |
depsjdy = 0.0d0 |
823 |
< |
depsjdz = 0.0d0 |
824 |
< |
depsjdux = 0.0d0 |
825 |
< |
depsjduy = 0.0d0 |
826 |
< |
depsjduz = 0.0d0 |
809 |
> |
dsigmajdx = 0.0_dp |
810 |
> |
dsigmajdy = 0.0_dp |
811 |
> |
dsigmajdz = 0.0_dp |
812 |
> |
dsigmajdux = 0.0_dp |
813 |
> |
dsigmajduy = 0.0_dp |
814 |
> |
dsigmajduz = 0.0_dp |
815 |
> |
dsjdx = 0.0_dp |
816 |
> |
dsjdy = 0.0_dp |
817 |
> |
dsjdz = 0.0_dp |
818 |
> |
dsjdux = 0.0_dp |
819 |
> |
dsjduy = 0.0_dp |
820 |
> |
dsjduz = 0.0_dp |
821 |
> |
depsjdx = 0.0_dp |
822 |
> |
depsjdy = 0.0_dp |
823 |
> |
depsjdz = 0.0_dp |
824 |
> |
depsjdux = 0.0_dp |
825 |
> |
depsjduy = 0.0_dp |
826 |
> |
depsjduz = 0.0_dp |
827 |
|
else |
828 |
|
|
829 |
|
#ifdef IS_MPI |
857 |
|
|
858 |
|
dctjdx = - zj * xj / r3 |
859 |
|
dctjdy = - zj * yj / r3 |
860 |
< |
dctjdz = 1.0d0 / rij - zj2 / r3 |
860 |
> |
dctjdz = 1.0_dp / rij - zj2 / r3 |
861 |
|
dctjdux = yj / rij !- (zi * xj2) / r3 |
862 |
|
dctjduy = -xj / rij !- (zj * yj2) / r3 |
863 |
< |
dctjduz = 0.0d0 !zj / rij - (zj2 * zj) / r3 |
863 |
> |
dctjduz = 0.0_dp !zj / rij - (zj2 * zj) / r3 |
864 |
|
|
865 |
|
! this is an attempt to try to truncate the singularity when |
866 |
|
! sin(theta) is near 0.0: |
867 |
|
|
868 |
|
stj2 = 1.0_dp - ctj*ctj |
869 |
< |
if (dabs(stj2) .lt. 1.0d-12) then |
870 |
< |
projj = sqrt(rij * 1.0d-12) |
871 |
< |
dcpjdx = 1.0d0 / projj |
872 |
< |
dcpjdy = 0.0d0 |
869 |
> |
if (abs(stj2) .lt. 1.0e-12_dp) then |
870 |
> |
projj = sqrt(rij * 1.0e-12_dp) |
871 |
> |
dcpjdx = 1.0_dp / projj |
872 |
> |
dcpjdy = 0.0_dp |
873 |
|
dcpjdux = xj / projj |
874 |
< |
dcpjduy = 0.0d0 |
875 |
< |
dspjdx = 0.0d0 |
876 |
< |
dspjdy = 1.0d0 / projj |
877 |
< |
dspjdux = 0.0d0 |
874 |
> |
dcpjduy = 0.0_dp |
875 |
> |
dspjdx = 0.0_dp |
876 |
> |
dspjdy = 1.0_dp / projj |
877 |
> |
dspjdux = 0.0_dp |
878 |
|
dspjduy = yj / projj |
879 |
|
else |
880 |
|
projj = sqrt(xj2 + yj2) |
881 |
|
projj3 = projj*projj*projj |
882 |
< |
dcpjdx = 1.0d0 / projj - xj2 / projj3 |
882 |
> |
dcpjdx = 1.0_dp / projj - xj2 / projj3 |
883 |
|
dcpjdy = - xj * yj / projj3 |
884 |
|
dcpjdux = xj / projj - (xj2 * xj) / projj3 |
885 |
|
dcpjduy = - (xj * yj2) / projj3 |
886 |
|
dspjdx = - xj * yj / projj3 |
887 |
< |
dspjdy = 1.0d0 / projj - yj2 / projj3 |
887 |
> |
dspjdy = 1.0_dp / projj - yj2 / projj3 |
888 |
|
dspjdux = - (yj * xj2) / projj3 |
889 |
|
dspjduy = yj / projj - (yj2 * yj) / projj3 |
890 |
|
endif |
891 |
|
|
892 |
|
cpj = xj / projj |
893 |
< |
dcpjdz = 0.0d0 |
894 |
< |
dcpjduz = 0.0d0 |
893 |
> |
dcpjdz = 0.0_dp |
894 |
> |
dcpjduz = 0.0_dp |
895 |
|
|
896 |
|
spj = yj / projj |
897 |
< |
dspjdz = 0.0d0 |
898 |
< |
dspjduz = 0.0d0 |
897 |
> |
dspjdz = 0.0_dp |
898 |
> |
dspjduz = 0.0_dp |
899 |
|
|
900 |
|
|
901 |
|
! write(*,*) 'dcpdu = ' ,dcpidux, dcpiduy, dcpiduz |
909 |
|
call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, MMAX, & |
910 |
|
CHEBYSHEV_UN, um_j, dum_j) |
911 |
|
|
912 |
< |
sigma_j = 0.0d0 |
913 |
< |
s_j = 0.0d0 |
914 |
< |
eps_j = 0.0d0 |
915 |
< |
dsigmajdx = 0.0d0 |
916 |
< |
dsigmajdy = 0.0d0 |
917 |
< |
dsigmajdz = 0.0d0 |
918 |
< |
dsigmajdux = 0.0d0 |
919 |
< |
dsigmajduy = 0.0d0 |
920 |
< |
dsigmajduz = 0.0d0 |
921 |
< |
dsjdx = 0.0d0 |
922 |
< |
dsjdy = 0.0d0 |
923 |
< |
dsjdz = 0.0d0 |
924 |
< |
dsjdux = 0.0d0 |
925 |
< |
dsjduy = 0.0d0 |
926 |
< |
dsjduz = 0.0d0 |
927 |
< |
depsjdx = 0.0d0 |
928 |
< |
depsjdy = 0.0d0 |
929 |
< |
depsjdz = 0.0d0 |
930 |
< |
depsjdux = 0.0d0 |
931 |
< |
depsjduy = 0.0d0 |
932 |
< |
depsjduz = 0.0d0 |
912 |
> |
sigma_j = 0.0_dp |
913 |
> |
s_j = 0.0_dp |
914 |
> |
eps_j = 0.0_dp |
915 |
> |
dsigmajdx = 0.0_dp |
916 |
> |
dsigmajdy = 0.0_dp |
917 |
> |
dsigmajdz = 0.0_dp |
918 |
> |
dsigmajdux = 0.0_dp |
919 |
> |
dsigmajduy = 0.0_dp |
920 |
> |
dsigmajduz = 0.0_dp |
921 |
> |
dsjdx = 0.0_dp |
922 |
> |
dsjdy = 0.0_dp |
923 |
> |
dsjdz = 0.0_dp |
924 |
> |
dsjdux = 0.0_dp |
925 |
> |
dsjduy = 0.0_dp |
926 |
> |
dsjduz = 0.0_dp |
927 |
> |
depsjdx = 0.0_dp |
928 |
> |
depsjdy = 0.0_dp |
929 |
> |
depsjdz = 0.0_dp |
930 |
> |
depsjdux = 0.0_dp |
931 |
> |
depsjduy = 0.0_dp |
932 |
> |
depsjduz = 0.0_dp |
933 |
|
|
934 |
|
do lm = 1, ShapeMap%Shapes(st2)%nContactFuncs |
935 |
|
l = ShapeMap%Shapes(st2)%ContactFuncLValue(lm) |
1099 |
|
!!$ write(*,*) 'dsidu = ', dsidux, dsiduy, dsiduz |
1100 |
|
!!$ write(*,*) 'dsigidu = ', dsigmaidux, dsigmaiduy, dsigmaiduz |
1101 |
|
!!$ write(*,*) sigma_j, ' is sigma j; ', s_j, ' is s j; ', eps_j, ' is eps j' |
1102 |
< |
depsdxi = eps_j * depsidx / (2.0d0 * eps) |
1103 |
< |
depsdyi = eps_j * depsidy / (2.0d0 * eps) |
1104 |
< |
depsdzi = eps_j * depsidz / (2.0d0 * eps) |
1105 |
< |
depsduxi = eps_j * depsidux / (2.0d0 * eps) |
1106 |
< |
depsduyi = eps_j * depsiduy / (2.0d0 * eps) |
1107 |
< |
depsduzi = eps_j * depsiduz / (2.0d0 * eps) |
1102 |
> |
depsdxi = eps_j * depsidx / (2.0_dp * eps) |
1103 |
> |
depsdyi = eps_j * depsidy / (2.0_dp * eps) |
1104 |
> |
depsdzi = eps_j * depsidz / (2.0_dp * eps) |
1105 |
> |
depsduxi = eps_j * depsidux / (2.0_dp * eps) |
1106 |
> |
depsduyi = eps_j * depsiduy / (2.0_dp * eps) |
1107 |
> |
depsduzi = eps_j * depsiduz / (2.0_dp * eps) |
1108 |
|
|
1109 |
< |
depsdxj = eps_i * depsjdx / (2.0d0 * eps) |
1110 |
< |
depsdyj = eps_i * depsjdy / (2.0d0 * eps) |
1111 |
< |
depsdzj = eps_i * depsjdz / (2.0d0 * eps) |
1112 |
< |
depsduxj = eps_i * depsjdux / (2.0d0 * eps) |
1113 |
< |
depsduyj = eps_i * depsjduy / (2.0d0 * eps) |
1114 |
< |
depsduzj = eps_i * depsjduz / (2.0d0 * eps) |
1109 |
> |
depsdxj = eps_i * depsjdx / (2.0_dp * eps) |
1110 |
> |
depsdyj = eps_i * depsjdy / (2.0_dp * eps) |
1111 |
> |
depsdzj = eps_i * depsjdz / (2.0_dp * eps) |
1112 |
> |
depsduxj = eps_i * depsjdux / (2.0_dp * eps) |
1113 |
> |
depsduyj = eps_i * depsjduy / (2.0_dp * eps) |
1114 |
> |
depsduzj = eps_i * depsjduz / (2.0_dp * eps) |
1115 |
|
|
1116 |
|
!!$ write(*,*) 'depsidu = ', depsidux, depsiduy, depsiduz |
1117 |
|
|
1147 |
|
rt12 = rt6*rt6 |
1148 |
|
rt126 = rt12 - rt6 |
1149 |
|
|
1150 |
< |
pot_temp = 4.0d0 * eps * rt126 |
1150 |
> |
pot_temp = 4.0_dp * eps * rt126 |
1151 |
|
|
1152 |
|
vpair = vpair + pot_temp |
1153 |
|
if (do_pot) then |
1154 |
|
#ifdef IS_MPI |
1155 |
< |
pot_row(atom1) = pot_row(atom1) + 0.5d0*pot_temp*sw |
1156 |
< |
pot_col(atom2) = pot_col(atom2) + 0.5d0*pot_temp*sw |
1155 |
> |
pot_row(VDW_POT,atom1) = pot_row(VDW_POT,atom1) + 0.5_dp*pot_temp*sw |
1156 |
> |
pot_col(VDW_POT,atom2) = pot_col(VDW_POT,atom2) + 0.5_dp*pot_temp*sw |
1157 |
|
#else |
1158 |
|
pot = pot + pot_temp*sw |
1159 |
|
#endif |
1161 |
|
|
1162 |
|
!!$ write(*,*) 'drtdu, depsdu = ', drtduxi, depsduxi |
1163 |
|
|
1164 |
< |
dvdxi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdxi + 4.0d0*depsdxi*rt126 |
1165 |
< |
dvdyi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdyi + 4.0d0*depsdyi*rt126 |
1166 |
< |
dvdzi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdzi + 4.0d0*depsdzi*rt126 |
1167 |
< |
dvduxi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduxi + 4.0d0*depsduxi*rt126 |
1168 |
< |
dvduyi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduyi + 4.0d0*depsduyi*rt126 |
1169 |
< |
dvduzi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduzi + 4.0d0*depsduzi*rt126 |
1164 |
> |
dvdxi = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtdxi + 4.0_dp*depsdxi*rt126 |
1165 |
> |
dvdyi = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtdyi + 4.0_dp*depsdyi*rt126 |
1166 |
> |
dvdzi = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtdzi + 4.0_dp*depsdzi*rt126 |
1167 |
> |
dvduxi = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtduxi + 4.0_dp*depsduxi*rt126 |
1168 |
> |
dvduyi = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtduyi + 4.0_dp*depsduyi*rt126 |
1169 |
> |
dvduzi = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtduzi + 4.0_dp*depsduzi*rt126 |
1170 |
|
|
1171 |
< |
dvdxj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdxj + 4.0d0*depsdxj*rt126 |
1172 |
< |
dvdyj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdyj + 4.0d0*depsdyj*rt126 |
1173 |
< |
dvdzj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdzj + 4.0d0*depsdzj*rt126 |
1174 |
< |
dvduxj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduxj + 4.0d0*depsduxj*rt126 |
1175 |
< |
dvduyj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduyj + 4.0d0*depsduyj*rt126 |
1176 |
< |
dvduzj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduzj + 4.0d0*depsduzj*rt126 |
1171 |
> |
dvdxj = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtdxj + 4.0_dp*depsdxj*rt126 |
1172 |
> |
dvdyj = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtdyj + 4.0_dp*depsdyj*rt126 |
1173 |
> |
dvdzj = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtdzj + 4.0_dp*depsdzj*rt126 |
1174 |
> |
dvduxj = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtduxj + 4.0_dp*depsduxj*rt126 |
1175 |
> |
dvduyj = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtduyj + 4.0_dp*depsduyj*rt126 |
1176 |
> |
dvduzj = 24.0_dp*eps*(2.0_dp*rt11 - rt5)*drtduzj + 4.0_dp*depsduzj*rt126 |
1177 |
|
!!$ write(*,*) 'drtduxi = ', drtduxi, ' depsduxi = ', depsduxi |
1178 |
|
! do the torques first since they are easy: |
1179 |
|
! remember that these are still in the body fixed axes |
1180 |
|
|
1181 |
< |
txi = 0.0d0 |
1182 |
< |
tyi = 0.0d0 |
1183 |
< |
tzi = 0.0d0 |
1181 |
> |
txi = 0.0_dp |
1182 |
> |
tyi = 0.0_dp |
1183 |
> |
tzi = 0.0_dp |
1184 |
|
|
1185 |
< |
txj = 0.0d0 |
1186 |
< |
tyj = 0.0d0 |
1187 |
< |
tzj = 0.0d0 |
1185 |
> |
txj = 0.0_dp |
1186 |
> |
tyj = 0.0_dp |
1187 |
> |
tzj = 0.0_dp |
1188 |
|
|
1189 |
|
txi = (dvduyi - dvduzi) * sw |
1190 |
|
tyi = (dvduzi - dvduxi) * sw |
1343 |
|
end DO |
1344 |
|
|
1345 |
|
! start with 0,0: |
1346 |
< |
PLM(0,0) = 1.0D0 |
1346 |
> |
PLM(0,0) = 1.0_DP |
1347 |
|
|
1348 |
|
! x = +/- 1 functions are easy: |
1349 |
< |
IF (abs(X).EQ.1.0D0) THEN |
1349 |
> |
IF (abs(X).EQ.1.0_DP) THEN |
1350 |
|
DO I = 1, m |
1351 |
|
PLM(0, I) = X**I |
1352 |
< |
DLM(0, I) = 0.5D0*I*(I+1.0D0)*X**(I+1) |
1352 |
> |
DLM(0, I) = 0.5_DP*I*(I+1.0_DP)*X**(I+1) |
1353 |
|
end DO |
1354 |
|
DO J = 1, m |
1355 |
|
DO I = 1, l |
1356 |
|
IF (I.EQ.1) THEN |
1357 |
|
DLM(I, J) = 1.0D+300 |
1358 |
|
ELSE IF (I.EQ.2) THEN |
1359 |
< |
DLM(I, J) = -0.25D0*(J+2)*(J+1)*J*(J-1)*X**(J+1) |
1359 |
> |
DLM(I, J) = -0.25_DP*(J+2)*(J+1)*J*(J-1)*X**(J+1) |
1360 |
|
ENDIF |
1361 |
|
end DO |
1362 |
|
end DO |
1364 |
|
ENDIF |
1365 |
|
|
1366 |
|
LS = 1 |
1367 |
< |
IF (abs(X).GT.1.0D0) LS = -1 |
1368 |
< |
XQ = sqrt(LS*(1.0D0-X*X)) |
1369 |
< |
XS = LS*(1.0D0-X*X) |
1367 |
> |
IF (abs(X).GT.1.0_DP) LS = -1 |
1368 |
> |
XQ = sqrt(LS*(1.0_DP-X*X)) |
1369 |
> |
XS = LS*(1.0_DP-X*X) |
1370 |
|
|
1371 |
|
DO I = 1, l |
1372 |
< |
PLM(I, I) = -LS*(2.0D0*I-1.0D0)*XQ*PLM(I-1, I-1) |
1372 |
> |
PLM(I, I) = -LS*(2.0_DP*I-1.0_DP)*XQ*PLM(I-1, I-1) |
1373 |
|
enddo |
1374 |
|
|
1375 |
|
DO I = 0, l |
1376 |
< |
PLM(I, I+1)=(2.0D0*I+1.0D0)*X*PLM(I, I) |
1376 |
> |
PLM(I, I+1)=(2.0_DP*I+1.0_DP)*X*PLM(I, I) |
1377 |
|
enddo |
1378 |
|
|
1379 |
|
DO I = 0, l |
1380 |
|
DO J = I+2, m |
1381 |
< |
PLM(I, J)=((2.0D0*J-1.0D0)*X*PLM(I,J-1) - & |
1382 |
< |
(I+J-1.0D0)*PLM(I,J-2))/(J-I) |
1381 |
> |
PLM(I, J)=((2.0_DP*J-1.0_DP)*X*PLM(I,J-1) - & |
1382 |
> |
(I+J-1.0_DP)*PLM(I,J-2))/(J-I) |
1383 |
|
end DO |
1384 |
|
end DO |
1385 |
|
|
1386 |
< |
DLM(0, 0)=0.0D0 |
1386 |
> |
DLM(0, 0)=0.0_DP |
1387 |
|
DO J = 1, m |
1388 |
|
DLM(0, J)=LS*J*(PLM(0,J-1)-X*PLM(0,J))/XS |
1389 |
|
end DO |
1390 |
|
|
1391 |
|
DO I = 1, l |
1392 |
|
DO J = I, m |
1393 |
< |
DLM(I,J) = LS*I*X*PLM(I, J)/XS + (J+I)*(J-I+1.0D0)/XQ*PLM(I-1, J) |
1393 |
> |
DLM(I,J) = LS*I*X*PLM(I, J)/XS + (J+I)*(J-I+1.0_DP)/XQ*PLM(I-1, J) |
1394 |
|
end DO |
1395 |
|
end DO |
1396 |
|
|
1419 |
|
! The original Fortran77 codes can be found here: |
1420 |
|
! http://iris-lee3.ece.uiuc.edu/~jjin/routines/routines.html |
1421 |
|
|
1422 |
< |
real(kind=8), intent(in) :: x |
1422 |
> |
real(kind=dp), intent(in) :: x |
1423 |
|
integer, intent(in):: m, mmax |
1424 |
|
integer, intent(in):: function_type |
1425 |
< |
real(kind=8), dimension(0:mmax), intent(inout) :: pl, dpl |
1425 |
> |
real(kind=dp), dimension(0:mmax), intent(inout) :: pl, dpl |
1426 |
|
|
1427 |
< |
real(kind=8) :: a, b, c, y0, y1, dy0, dy1, yn, dyn |
1427 |
> |
real(kind=dp) :: a, b, c, y0, y1, dy0, dy1, yn, dyn |
1428 |
|
integer :: k |
1429 |
|
|
1430 |
< |
A = 2.0D0 |
1431 |
< |
B = 0.0D0 |
1432 |
< |
C = 1.0D0 |
1433 |
< |
Y0 = 1.0D0 |
1434 |
< |
Y1 = 2.0D0*X |
1435 |
< |
DY0 = 0.0D0 |
1436 |
< |
DY1 = 2.0D0 |
1437 |
< |
PL(0) = 1.0D0 |
1438 |
< |
PL(1) = 2.0D0*X |
1439 |
< |
DPL(0) = 0.0D0 |
1440 |
< |
DPL(1) = 2.0D0 |
1430 |
> |
A = 2.0_DP |
1431 |
> |
B = 0.0_DP |
1432 |
> |
C = 1.0_DP |
1433 |
> |
Y0 = 1.0_DP |
1434 |
> |
Y1 = 2.0_DP*X |
1435 |
> |
DY0 = 0.0_DP |
1436 |
> |
DY1 = 2.0_DP |
1437 |
> |
PL(0) = 1.0_DP |
1438 |
> |
PL(1) = 2.0_DP*X |
1439 |
> |
DPL(0) = 0.0_DP |
1440 |
> |
DPL(1) = 2.0_DP |
1441 |
|
IF (function_type.EQ.CHEBYSHEV_TN) THEN |
1442 |
|
Y1 = X |
1443 |
< |
DY1 = 1.0D0 |
1443 |
> |
DY1 = 1.0_DP |
1444 |
|
PL(1) = X |
1445 |
< |
DPL(1) = 1.0D0 |
1445 |
> |
DPL(1) = 1.0_DP |
1446 |
|
ELSE IF (function_type.EQ.LAGUERRE) THEN |
1447 |
< |
Y1 = 1.0D0-X |
1448 |
< |
DY1 = -1.0D0 |
1449 |
< |
PL(1) = 1.0D0-X |
1450 |
< |
DPL(1) = -1.0D0 |
1447 |
> |
Y1 = 1.0_DP-X |
1448 |
> |
DY1 = -1.0_DP |
1449 |
> |
PL(1) = 1.0_DP-X |
1450 |
> |
DPL(1) = -1.0_DP |
1451 |
|
ENDIF |
1452 |
|
DO K = 2, m |
1453 |
|
IF (function_type.EQ.LAGUERRE) THEN |
1454 |
< |
A = -1.0D0/K |
1455 |
< |
B = 2.0D0+A |
1456 |
< |
C = 1.0D0+A |
1454 |
> |
A = -1.0_DP/K |
1455 |
> |
B = 2.0_DP+A |
1456 |
> |
C = 1.0_DP+A |
1457 |
|
ELSE IF (function_type.EQ.HERMITE) THEN |
1458 |
< |
C = 2.0D0*(K-1.0D0) |
1458 |
> |
C = 2.0_DP*(K-1.0_DP) |
1459 |
|
ENDIF |
1460 |
|
YN = (A*X+B)*Y1-C*Y0 |
1461 |
|
DYN = A*Y1+(A*X+B)*DY1-C*DY0 |