TABLE OF CONTENTS
- ABINIT/align_u_matrices
- ABINIT/calc_b_matrix
- ABINIT/calc_btinv_matrix
- ABINIT/dang_d1
- ABINIT/dang_d2
- ABINIT/dbond_length_d1
- ABINIT/ddihedral_d1
- ABINIT/ddihedral_d2
- ABINIT/deloc2xcart
- ABINIT/gred2gdeloc
- ABINIT/m_pred_delocint
- ABINIT/pred_delocint
- ABINIT/xcart2deloc
- ABINIT/xfh_recover_deloc
ABINIT/align_u_matrices [ Functions ]
NAME
align_u_matrices
FUNCTION
INPUTS
OUTPUT
SOURCE
1570 subroutine align_u_matrices(natom,ninternal,u_matrix,u_matrix_old,s_matrix,f_eigs) 1571 1572 !Arguments ------------------------------------ 1573 !scalars 1574 integer,intent(in) :: ninternal,natom 1575 !arrays 1576 real(dp),intent(in) :: u_matrix_old(ninternal,3*(natom-1)) 1577 real(dp),intent(inout) :: f_eigs(3*natom) 1578 real(dp),intent(inout) :: s_matrix(3*natom,3*natom) 1579 real(dp),intent(inout) :: u_matrix(ninternal,3*(natom-1)) 1580 1581 !Local variables ------------------------------ 1582 !scalars 1583 integer :: ii,iint1,imax 1584 real(dp) :: ss 1585 !arrays 1586 integer :: eigv_flag(3*(natom-1)),eigv_ind(3*(natom-1)) 1587 real(dp) :: tmps(3*natom,3*natom) 1588 real(dp) :: tmpu(ninternal,3*(natom-1)) 1589 real(dp) :: tmpf(3*natom) 1590 1591 !****************************************************************** 1592 1593 eigv_flag(:) = 0 1594 eigv_ind(:) = 0 1595 1596 !just permit a change in sign 1597 do iint1=1,3*(natom-1) 1598 ss = zero 1599 do ii=1,ninternal 1600 ss = ss + u_matrix_old(ii,iint1)*u_matrix(ii,iint1) 1601 end do 1602 if (ss < -tol12) then 1603 imax = -iint1 1604 else 1605 imax = iint1 1606 end if 1607 eigv_ind(iint1) = imax 1608 eigv_flag(abs(imax)) = 1 1609 end do 1610 1611 tmpu(:,:) = u_matrix 1612 tmps(:,:) = s_matrix 1613 tmpf(:) = f_eigs 1614 !exchange eigenvectors... 1615 do iint1=1,3*(natom-1) 1616 ss = one 1617 if (eigv_ind(iint1) < 0) ss = -one 1618 1619 imax = abs(eigv_ind(iint1)) 1620 1621 tmpu(:,imax) = ss*u_matrix(:,iint1) 1622 1623 tmps(:,imax+3) = ss*s_matrix(:,iint1+3) 1624 1625 tmpf(imax+3) = f_eigs(iint1+3) 1626 end do 1627 1628 u_matrix(:,:) = tmpu(:,:) 1629 s_matrix(:,:) = tmps(:,:) 1630 f_eigs(:) = tmpf(:) 1631 1632 end subroutine align_u_matrices
ABINIT/calc_b_matrix [ Functions ]
NAME
calc_b_matrix
FUNCTION
calculate values of derivatives of internal coordinates as a function of cartesian ones = B matrix
INPUTS
angs= number of angles
bonds(2,2,nbond)=for a bond between iatom and jatom
bonds(1,1,nbond) = iatom
bonds(2,1,nbond) = icenter
bonds(1,2,nbond) = jatom
bonds(2,2,nbond) = irshift
carts(2,ncart)= index of total primitive internal, and atom (carts(2,:))
dihedrals(2,4,ndihed)=indexes to characterize dihedrals
nang(2,3,nang)=indexes to characterize angles
nbond=number of bonds
ncart=number of auxiliary cartesian atom coordinates (used for constraints)
ndihed= number of dihedrals
ninternal=nbond+nang+ndihed+ncart: number of internal coordinates
nrshift= dimension of rshift
rprimd(3,3)=dimensional real space primitive translations (bohr)
rshift(3,nrshift)=shift in xred that must be done to find all neighbors of
a given atom within a given number of neighboring shells
xcart(3,natom)=cartesian coordinates of atoms (bohr)
OUTPUT
b_matrix(ninternal,3*natom)=matrix of derivatives of internal coordinates wrt cartesians
SOURCE
887 subroutine calc_b_matrix(deloc,natom,rprimd,xcart,b_matrix) 888 889 !Arguments ------------------------------------ 890 !scalars 891 integer,intent(in) :: natom 892 type(delocint),intent(in) :: deloc 893 894 !arrays 895 real(dp),intent(in) :: rprimd(3,3),xcart(3,natom) 896 real(dp),intent(out) :: b_matrix(deloc%ninternal,3*natom) 897 898 !Local variables------------------------------- 899 !scalars 900 integer :: i1,i2,i3,i4,iang,ibond,icart,idihed,iprim,s1,s2,s3,s4 901 !arrays 902 real(dp) :: bb(3),r1(3),r2(3),r3(3),r4(3) 903 904 ! ************************************************************************* 905 906 iprim=0 907 b_matrix(:,:) = zero 908 909 do ibond=1,deloc%nbond 910 i1 = deloc%bonds(1,1,ibond) 911 s1 = deloc%bonds(2,1,ibond) 912 r1(:) = xcart(:,i1)+deloc%rshift(1,s1)*rprimd(:,1)& 913 & +deloc%rshift(2,s1)*rprimd(:,2)& 914 & +deloc%rshift(3,s1)*rprimd(:,3) 915 i2 = deloc%bonds(1,2,ibond) 916 s2 = deloc%bonds(2,2,ibond) 917 r2(:) = xcart(:,i2)+deloc%rshift(1,s2)*rprimd(:,1)& 918 & +deloc%rshift(2,s2)*rprimd(:,2)& 919 & +deloc%rshift(3,s2)*rprimd(:,3) 920 iprim=iprim+1 921 call dbond_length_d1(r1,r2,bb) 922 b_matrix(iprim,3*(i1-1)+1:3*i1) = b_matrix(iprim,3*(i1-1)+1:3*i1) + bb(:) 923 call dbond_length_d1(r2,r1,bb) 924 b_matrix(iprim,3*(i2-1)+1:3*i2) = b_matrix(iprim,3*(i2-1)+1:3*i2) + bb(:) 925 end do 926 927 !second: angle values (ang) 928 do iang=1,deloc%nang 929 i1 = deloc%angs(1,1,iang) 930 s1 = deloc%angs(2,1,iang) 931 r1(:) = xcart(:,i1)+deloc%rshift(1,s1)*rprimd(:,1)& 932 & +deloc%rshift(2,s1)*rprimd(:,2)& 933 & +deloc%rshift(3,s1)*rprimd(:,3) 934 i2 = deloc%angs(1,2,iang) 935 s2 = deloc%angs(2,2,iang) 936 r2(:) = xcart(:,i2)+deloc%rshift(1,s2)*rprimd(:,1)& 937 & +deloc%rshift(2,s2)*rprimd(:,2)& 938 & +deloc%rshift(3,s2)*rprimd(:,3) 939 i3 = deloc%angs(1,3,iang) 940 s3 = deloc%angs(2,3,iang) 941 r3(:) = xcart(:,i3)+deloc%rshift(1,s3)*rprimd(:,1)& 942 & +deloc%rshift(2,s3)*rprimd(:,2)& 943 & +deloc%rshift(3,s3)*rprimd(:,3) 944 iprim=iprim+1 945 call dang_d1(r1,r2,r3,bb) 946 b_matrix(iprim,3*(i1-1)+1:3*i1) = b_matrix(iprim,3*(i1-1)+1:3*i1) + bb(:) 947 call dang_d2(r1,r2,r3,bb) 948 b_matrix(iprim,3*(i2-1)+1:3*i2) = b_matrix(iprim,3*(i2-1)+1:3*i2) + bb(:) 949 call dang_d1(r3,r2,r1,bb) 950 b_matrix(iprim,3*(i3-1)+1:3*i3) = b_matrix(iprim,3*(i3-1)+1:3*i3) + bb(:) 951 end do 952 953 !third: dihedral values 954 do idihed=1,deloc%ndihed 955 i1 = deloc%dihedrals(1,1,idihed) 956 s1 = deloc%dihedrals(2,1,idihed) 957 r1(:) = xcart(:,i1)+deloc%rshift(1,s1)*rprimd(:,1)& 958 & +deloc%rshift(2,s1)*rprimd(:,2)& 959 & +deloc%rshift(3,s1)*rprimd(:,3) 960 i2 = deloc%dihedrals(1,2,idihed) 961 s2 = deloc%dihedrals(2,2,idihed) 962 r2(:) = xcart(:,i2)+deloc%rshift(1,s2)*rprimd(:,1)& 963 & +deloc%rshift(2,s2)*rprimd(:,2)& 964 & +deloc%rshift(3,s2)*rprimd(:,3) 965 i3 = deloc%dihedrals(1,3,idihed) 966 s3 = deloc%dihedrals(2,3,idihed) 967 r3(:) = xcart(:,i3)+deloc%rshift(1,s3)*rprimd(:,1)& 968 & +deloc%rshift(2,s3)*rprimd(:,2)& 969 & +deloc%rshift(3,s3)*rprimd(:,3) 970 i4 = deloc%dihedrals(1,4,idihed) 971 s4 = deloc%dihedrals(2,4,idihed) 972 r4(:) = xcart(:,i4)+deloc%rshift(1,s4)*rprimd(:,1)& 973 & +deloc%rshift(2,s4)*rprimd(:,2)& 974 & +deloc%rshift(3,s4)*rprimd(:,3) 975 ! write(std_out,*) 'dihed ',idihed 976 ! write(std_out,*) r1 977 ! write(std_out,*) r2 978 ! write(std_out,*) r3 979 ! write(std_out,*) r4 980 981 iprim=iprim+1 982 call ddihedral_d1(r1,r2,r3,r4,bb) 983 b_matrix(iprim,3*(i1-1)+1:3*i1) = b_matrix(iprim,3*(i1-1)+1:3*i1) + bb(:) 984 call ddihedral_d2(r1,r2,r3,r4,bb) 985 b_matrix(iprim,3*(i2-1)+1:3*i2) = b_matrix(iprim,3*(i2-1)+1:3*i2) + bb(:) 986 call ddihedral_d2(r4,r3,r2,r1,bb) 987 b_matrix(iprim,3*(i3-1)+1:3*i3) = b_matrix(iprim,3*(i3-1)+1:3*i3) + bb(:) 988 call ddihedral_d1(r4,r3,r2,r1,bb) 989 b_matrix(iprim,3*(i4-1)+1:3*i4) = b_matrix(iprim,3*(i4-1)+1:3*i4) + bb(:) 990 end do 991 992 do icart=1,deloc%ncart 993 iprim=iprim+1 994 b_matrix(iprim,3*(deloc%carts(2,icart)-1)+deloc%carts(1,icart)) = & 995 & b_matrix(iprim,3*(deloc%carts(2,icart)-1)+deloc%carts(1,icart)) + one 996 end do 997 998 !DEBUG 999 ! write (200,*) 'calc_b_matrix : b_matrix = ' 1000 ! do iprim=1,deloc%ninternal 1001 ! do i1=1, 3*natom 1002 ! write (200,'(E16.6,2x)',ADVANCE='NO') b_matrix(iprim,i1) 1003 ! end do 1004 ! write (200,*) 1005 ! end do 1006 !ENDDEBUG 1007 1008 end subroutine calc_b_matrix
ABINIT/calc_btinv_matrix [ Functions ]
NAME
calc_btinv_matrix
FUNCTION
INPUTS
OUTPUT
NOTES
bt_inv_matrix is inverse transpose of the delocalized
coordinate B matrix. b_matrix is the primitive internal B matrix
SOURCE
1481 subroutine calc_btinv_matrix(b_matrix,natom,ninternal,bt_inv_matrix,u_matrix) 1482 1483 !Arguments ------------------------------------ 1484 integer,intent(in) :: ninternal,natom 1485 real(dp),intent(in) :: b_matrix(ninternal,3*natom) 1486 real(dp),intent(out) :: bt_inv_matrix(3*(natom-1),3*natom) 1487 real(dp),intent(inout) :: u_matrix(ninternal,3*(natom-1)) 1488 1489 !Local variables ------------------------------------ 1490 !scalars 1491 integer :: ii,info,lwork 1492 !arrays 1493 real(dp) :: f_eigs(3*natom),f_matrix(3*natom,3*natom) 1494 real(dp) :: s_matrix(3*natom,3*natom) 1495 real(dp) :: s_red(3*natom,3*(natom-1)) 1496 real(dp) :: u_matrix_old(ninternal,3*(natom-1)) 1497 real(dp),allocatable :: work(:) 1498 1499 !****************************************************************** 1500 1501 !f matrix = B^{T} B 1502 call dgemm('T','N',3*natom,3*natom,ninternal,one,& 1503 & b_matrix,ninternal,b_matrix,ninternal,zero,f_matrix,3*natom) 1504 1505 lwork = max(1,3*3*natom-1) 1506 ABI_MALLOC(work,(lwork)) 1507 s_matrix(:,:) = f_matrix(:,:) 1508 1509 call dsyev('V','L',3*natom,s_matrix,3*natom,f_eigs,work,lwork,info) 1510 1511 ABI_FREE(work) 1512 1513 if (abs(f_eigs(1)) + abs(f_eigs(2)) + abs(f_eigs(3)) > tol10 ) then 1514 write(std_out,*) 'Error: 3 lowest eigenvalues are not zero' 1515 write(std_out,*) ' internal coordinates do NOT span the full degrees of freedom !' 1516 write(std_out,'(6E16.6)') f_eigs 1517 ABI_ERROR("Aborting now") 1518 end if 1519 if ( abs(f_eigs(4)) < tol10 ) then 1520 write(std_out,*) 'Error: fourth eigenvalue is zero' 1521 write(std_out,*) ' internal coordinates do NOT span the full degrees of freedom !' 1522 write(std_out,'(6E16.6)') f_eigs 1523 ABI_ERROR("Aborting now") 1524 end if 1525 1526 !calculate U matrix from U = B * S_red * lambda^{-1/2} 1527 do ii=1,3*(natom-1) 1528 s_red(:,ii) = s_matrix(:,ii+3)/sqrt(f_eigs(ii+3)) 1529 end do 1530 1531 u_matrix_old(:,:) = u_matrix(:,:) 1532 1533 call dgemm('N','N',ninternal,3*(natom-1),3*natom,one,& 1534 & b_matrix,ninternal,s_red,3*natom,zero,u_matrix,ninternal) 1535 1536 1537 !align eigenvectors, to preserve a form of continuity in convergences 1538 !!!! eigenvalues are no longer in increasing order!!! but only s_red is reordered 1539 !so that btinv is correct. 1540 call align_u_matrices(natom,ninternal,u_matrix,u_matrix_old,s_matrix,f_eigs) 1541 1542 !calculate B_deloc^{-1} matrix for transformation of forces to deloc coord. 1543 !(B^{T}_deloc)^{-1} = (B_deloc B^{T}_deloc)^{-1} B_deloc = lambda^{-3/2} S^{T} F 1544 != ( S lambda^{3/2} )^{T} F 1545 1546 !! DEFINITION 1547 !! real(dp),intent(out) :: bt_inv_matrix(3*(natom-1),3*natom) 1548 1549 !even better: B_deloc^{-1} = lambda^{-1/2} S^{T} 1550 do ii=1,3*(natom-1) 1551 ! s_red(:,ii) = s_matrix(:,ii+3)*sqrt(f_eigs(ii+3)) 1552 bt_inv_matrix(ii,:) = s_matrix(:,ii+3)/sqrt(f_eigs(ii+3)) 1553 end do 1554 1555 end subroutine calc_btinv_matrix
ABINIT/dang_d1 [ Functions ]
NAME
dang_d1
FUNCTION
SOURCE
1048 subroutine dang_d1(r1,r2,r3,bb) 1049 1050 !Arguments ------------------------------------ 1051 !arrays 1052 real(dp),intent(in) :: r1(3),r2(3),r3(3) 1053 real(dp),intent(out) :: bb(3) 1054 1055 !Local variables ------------------------------ 1056 !scalars 1057 real(dp) :: cos_ang,n1,n1232,n2,tmp 1058 !arrays 1059 real(dp) :: cp1232(3),rpt(3),rpt12(3),rpt32(3) 1060 1061 !************************************************************************ 1062 n1=bond_length(r1,r2) 1063 n2=bond_length(r3,r2) 1064 1065 rpt12(:) = r1(:)-r2(:) 1066 rpt32(:) = r3(:)-r2(:) 1067 1068 cos_ang = (rpt12(1)*rpt32(1)+rpt12(2)*rpt32(2)+rpt12(3)*rpt32(3))/n1/n2 1069 if (cos_ang > one - epsilon(one)*two) then 1070 cos_ang = one 1071 else if(cos_ang < -one + epsilon(one)*two) then 1072 cos_ang = -one 1073 end if 1074 1075 rpt(:) = rpt32(:)/n1/n2 - rpt12(:)*cos_ang/n1/n1 1076 1077 tmp = sqrt(one-cos_ang**2) 1078 bb(:) = zero 1079 if (tmp > epsilon(one)) then 1080 bb(:) = rpt(:) * (-one)/tmp 1081 end if 1082 1083 !TEST: version from MOLECULAR VIBRATIONS EB Wilson 1084 call acrossb(rpt12,rpt32,cp1232) 1085 n1232 = sqrt(cp1232(1)**2+cp1232(2)**2+cp1232(3)**2) 1086 rpt(:) = (cos_ang*rpt12(:)*n2/n1 - rpt32(:))/n1232 1087 if (abs(bb(1)-rpt(1))+abs(bb(2)-rpt(2))+abs(bb(3)-rpt(3)) > tol10) then 1088 write(std_out,*) 'Compare bb ang 1 : ' 1089 write(std_out,*) bb(:), rpt(:), bb(:)-rpt(:) 1090 end if 1091 bb(:) = rpt(:) 1092 1093 end subroutine dang_d1
ABINIT/dang_d2 [ Functions ]
NAME
dang_d2
FUNCTION
SOURCE
1106 subroutine dang_d2(r1,r2,r3,bb) 1107 1108 !Arguments ------------------------------------ 1109 !arrays 1110 real(dp),intent(in) :: r1(3),r2(3),r3(3) 1111 real(dp),intent(out) :: bb(3) 1112 1113 !Local variables ------------------------------ 1114 !scalars 1115 real(dp) :: cos_ang,n1,n1232,n2,tmp 1116 !arrays 1117 real(dp) :: cp1232(3),rpt(3),rpt12(3),rpt32(3) 1118 1119 !************************************************************************ 1120 n1=bond_length(r1,r2) 1121 n2=bond_length(r3,r2) 1122 1123 rpt12(:) = r1(:)-r2(:) 1124 rpt32(:) = r3(:)-r2(:) 1125 1126 cos_ang = (rpt12(1)*rpt32(1)+rpt12(2)*rpt32(2)+rpt12(3)*rpt32(3))/n1/n2 1127 if (cos_ang > one - epsilon(one)*two) then 1128 cos_ang = one 1129 else if(cos_ang < -one + epsilon(one)*two) then 1130 cos_ang = -one 1131 end if 1132 1133 rpt(:) = -rpt32(:)/n1/n2 - rpt12(:)/n1/n2 & 1134 & + rpt12(:)*cos_ang/n1/n1 + rpt32(:)*cos_ang/n2/n2 1135 1136 tmp = sqrt(one-cos_ang**2) 1137 bb(:) = zero 1138 if (tmp > tol12) then 1139 bb(:) = rpt(:) * (-one)/tmp 1140 end if 1141 1142 !TEST: version from MOLECULAR VIBRATIONS EB Wilson 1143 call acrossb(rpt12,rpt32,cp1232) 1144 n1232 = sqrt(cp1232(1)**2+cp1232(2)**2+cp1232(3)**2) 1145 rpt(:) = ((n1-n2*cos_ang)*rpt12(:)/n1 + (n2-n1*cos_ang)*rpt32(:)/n2) / n1232 1146 if (abs(bb(1)-rpt(1))+abs(bb(2)-rpt(2))+abs(bb(3)-rpt(3)) > tol10) then 1147 write(std_out,*) 'Compare bb ang 2 : ' 1148 write(std_out,*) bb(:), rpt(:), bb(:)-rpt(:) 1149 end if 1150 bb(:) = rpt(:) 1151 1152 end subroutine dang_d2
ABINIT/dbond_length_d1 [ Functions ]
NAME
dbond_length_d1
FUNCTION
SOURCE
1020 subroutine dbond_length_d1(r1,r2,bb) 1021 1022 !Arguments ------------------------------------ 1023 !arrays 1024 real(dp),intent(in) :: r1(3),r2(3) 1025 real(dp),intent(out) :: bb(3) 1026 1027 !Local variables ------------------------------ 1028 !arrays 1029 real(dp) :: rpt(3) 1030 1031 !************************************************************************ 1032 rpt(:) = r1(:)-r2(:) 1033 bb(:) = rpt(:)/bond_length(r1,r2) 1034 1035 end subroutine dbond_length_d1
ABINIT/ddihedral_d1 [ Functions ]
NAME
ddihedral_d1
FUNCTION
SOURCE
1164 subroutine ddihedral_d1(r1,r2,r3,r4,bb) 1165 1166 !Arguments ------------------------------------ 1167 !arrays 1168 real(dp),intent(in) :: r1(3),r2(3),r3(3),r4(3) 1169 real(dp),intent(out) :: bb(3) 1170 1171 !Local variables ------------------------------------ 1172 !scalars 1173 real(dp) :: cos_dihedral,dih_sign,n1,n2,n23,sin_dihedral,tmp 1174 !arrays 1175 real(dp) :: cp1232(3),cp32_1232(3),cp32_3432(3),cp3432(3),cpcp(3),rpt(3) 1176 real(dp) :: rpt12(3),rpt32(3),rpt34(3) 1177 1178 !****************************************************************** 1179 rpt12(:) = r1(:)-r2(:) 1180 rpt32(:) = r3(:)-r2(:) 1181 rpt34(:) = r3(:)-r4(:) 1182 1183 call acrossb(rpt12,rpt32,cp1232) 1184 call acrossb(rpt34,rpt32,cp3432) 1185 1186 !DEBUG 1187 !write(std_out,*) ' cos_dihedral : cp1232 = ', cp1232 1188 !write(std_out,*) ' cos_dihedral : cp3432 = ', cp3432 1189 !ENDDEBUG 1190 1191 n1 = sqrt(cp1232(1)**2+cp1232(2)**2+cp1232(3)**2) 1192 n2 = sqrt(cp3432(1)**2+cp3432(2)**2+cp3432(3)**2) 1193 1194 cos_dihedral = (cp1232(1)*cp3432(1)+cp1232(2)*cp3432(2)+cp1232(3)*cp3432(3))/n1/n2 1195 if (cos_dihedral > one - epsilon(one)*two) then 1196 cos_dihedral = one 1197 else if(cos_dihedral < -one + epsilon(one)*two) then 1198 cos_dihedral = -one 1199 end if 1200 !we use complementary of standard angle, so 1201 !cos_dihedral = -cos_dihedral 1202 1203 call acrossb(cp1232,cp3432,cpcp) 1204 cpcp(:) = cpcp(:)/n1/n2 1205 !we use complementary of standard angle, but sin is invariant 1206 sin_dihedral = -(cpcp(1)*rpt32(1)+cpcp(2)*rpt32(2)+cpcp(3)*rpt32(3))& 1207 & /sqrt(rpt32(1)**2+rpt32(2)**2+rpt32(3)**2) 1208 dih_sign = one 1209 if (sin_dihedral < -epsilon(one)) then 1210 dih_sign = -one 1211 end if 1212 1213 !DEBUG 1214 !write(std_out,'(a,3E16.6)') 'ddihedral_d1 : cos abs(sin) dih_sign= ',& 1215 !& cos_dihedral,sin_dihedral,dih_sign 1216 !ENDDEBUG 1217 1218 !ddihedral_d1 = dih_sign* acos(cos_dihedral) 1219 call acrossb(rpt32,cp1232,cp32_1232) 1220 call acrossb(rpt32,cp3432,cp32_3432) 1221 1222 rpt(:) = cp32_3432(:)/n1/n2 - cp32_1232(:)/n1/n1 * cos_dihedral 1223 bb(:) = zero 1224 1225 !DEBUG 1226 !write(std_out,*) 'ddihedral_d1 cp1232 cp3432 = ',cp1232,cp3432,rpt32 1227 !write(std_out,*) 'ddihedral_d1 cp32_1232 cp32_3432 = ',cp32_1232,cp32_3432,cos_dihedral,n1,n2 1228 !write(std_out,*) 'ddihedral_d1 rpt = ',rpt 1229 !ENDDEBUG 1230 1231 tmp = sqrt(one-cos_dihedral**2) 1232 if (tmp > tol12) then 1233 ! we use complementary of standard angle, so cosine in acos has - sign, 1234 ! and it appears for the derivative 1235 bb(:) = -dih_sign * rpt(:) * (-one) / tmp 1236 else 1237 bb(:) = dih_sign * cp32_3432(:) / n1 / n2 / & 1238 & sqrt(cp32_3432(1)**2+cp32_3432(2)**2+cp32_3432(3)**2) 1239 end if 1240 1241 !TEST: version from MOLECULAR VIBRATIONS EB Wilson 1242 1243 n23 = sqrt(rpt32(1)*rpt32(1)+rpt32(2)*rpt32(2)+rpt32(3)*rpt32(3)) 1244 rpt(:) = cp1232(:)*n23/n1/n1 1245 !if (abs(bb(1)-rpt(1))+abs(bb(2)-rpt(2))+abs(bb(3)-rpt(3)) > tol10) then 1246 !write(std_out,*) 'Compare bb1 : ' 1247 !write(std_out,*) bb(:), rpt(:), bb(:)-rpt(:) 1248 !end if 1249 bb(:) = rpt(:) 1250 1251 end subroutine ddihedral_d1
ABINIT/ddihedral_d2 [ Functions ]
NAME
ddihedral_d2
FUNCTION
SOURCE
1263 subroutine ddihedral_d2(r1,r2,r3,r4,bb) 1264 1265 !Arguments ------------------------------------ 1266 !arrays 1267 real(dp),intent(in) :: r1(3),r2(3),r3(3),r4(3) 1268 real(dp),intent(out) :: bb(3) 1269 1270 !Local variables 1271 !scalars 1272 real(dp) :: cos_dihedral,dih_sign,n1,n2,n23,sin_dihedral,sp1232,sp3432,tmp 1273 !arrays 1274 real(dp) :: cp1232(3),cp1232_12(3),cp1232_34(3),cp32_1232(3),cp32_3432(3) 1275 real(dp) :: cp3432(3),cp3432_12(3),cp3432_34(3),cpcp(3),rpt(3),rpt12(3) 1276 real(dp) :: rpt32(3),rpt34(3) 1277 1278 ! ************************************************************************* 1279 rpt12(:) = r1(:)-r2(:) 1280 rpt32(:) = r3(:)-r2(:) 1281 rpt34(:) = r3(:)-r4(:) 1282 1283 call acrossb(rpt12,rpt32,cp1232) 1284 call acrossb(rpt34,rpt32,cp3432) 1285 1286 !DEBUG 1287 !write(std_out,*) ' cos_dihedral : cp1232 = ', cp1232 1288 !write(std_out,*) ' cos_dihedral : cp3432 = ', cp3432 1289 !ENDDEBUG 1290 1291 n1 = sqrt(cp1232(1)**2+cp1232(2)**2+cp1232(3)**2) 1292 n2 = sqrt(cp3432(1)**2+cp3432(2)**2+cp3432(3)**2) 1293 1294 cos_dihedral = (cp1232(1)*cp3432(1)+cp1232(2)*cp3432(2)+cp1232(3)*cp3432(3))/n1/n2 1295 if (cos_dihedral > one - epsilon(one)*two) then 1296 cos_dihedral = one 1297 else if(cos_dihedral < -one + epsilon(one)*two) then 1298 cos_dihedral = -one 1299 end if 1300 !we use complementary of standard angle, so 1301 !cos_dihedral = -cos_dihedral 1302 1303 call acrossb(cp1232,cp3432,cpcp) 1304 cpcp(:) = cpcp(:)/n1/n2 1305 !we use complementary of standard angle, but sin is invariant 1306 sin_dihedral = -(cpcp(1)*rpt32(1)+cpcp(2)*rpt32(2)+cpcp(3)*rpt32(3))& 1307 & /sqrt(rpt32(1)**2+rpt32(2)**2+rpt32(3)**2) 1308 dih_sign = one 1309 if (sin_dihedral < -tol12) then 1310 dih_sign = -one 1311 end if 1312 1313 !DEBUG 1314 !write(std_out,'(a,3E16.6)') 'ddihedral_d2 : cos abs(sin) dih_sign= ',& 1315 !& cos_dihedral,sin_dihedral,dih_sign 1316 !ENDDEBUG 1317 1318 !ddihedral_d2 = dih_sign* acos(cos_dihedral) 1319 call acrossb(rpt32,cp3432,cp32_3432) 1320 call acrossb(cp3432,rpt12,cp3432_12) 1321 call acrossb(cp1232,rpt34,cp1232_34) 1322 1323 call acrossb(rpt32,cp1232,cp32_1232) 1324 call acrossb(cp1232,rpt12,cp1232_12) 1325 call acrossb(cp3432,rpt34,cp3432_34) 1326 1327 rpt(:) = -(cp32_3432(:) + cp3432_12(:) + cp1232_34(:))/n1/n2 & 1328 & +cos_dihedral*(cp32_1232(:)/n1/n1 + cp1232_12(:)/n1/n1 + cp3432_34(:)/n2/n2) 1329 bb(:) = zero 1330 tmp = sqrt(one-cos_dihedral**2) 1331 if (tmp > tol12) then 1332 ! we use complementary of standard angle, so cosine in acos has - sign, 1333 ! and it appears for derivative 1334 bb(:) = -dih_sign * rpt(:) * (-one) / tmp 1335 else 1336 bb(:) = dih_sign * cos_dihedral * & 1337 & ( cp32_1232(:)/n1/n1/sqrt(cp32_1232(1)**2+cp32_1232(2)**2+cp32_1232(3)**2) & 1338 & +cp1232_12(:)/n1/n1/sqrt(cp1232_12(1)**2+cp1232_12(2)**2+cp1232_12(3)**2) & 1339 & +cp3432_34(:)/n2/n2/sqrt(cp3432_34(1)**2+cp3432_34(2)**2+cp3432_34(3)**2) ) 1340 end if 1341 1342 !TEST: version from MOLECULAR VIBRATIONS EB Wilson p. 61 1343 n23 = sqrt(rpt32(1)*rpt32(1)+rpt32(2)*rpt32(2)+rpt32(3)*rpt32(3)) 1344 sp1232 = rpt12(1)*rpt32(1)+rpt12(2)*rpt32(2)+rpt12(3)*rpt32(3) 1345 sp3432 = rpt34(1)*rpt32(1)+rpt34(2)*rpt32(2)+rpt34(3)*rpt32(3) 1346 1347 rpt(:) = -cp1232(:)*(n23-sp1232/n23)/n1/n1 - cp3432(:)*sp3432/n23/n2/n2 1348 !if (abs(bb(1)-rpt(1))+abs(bb(2)-rpt(2))+abs(bb(3)-rpt(3)) > tol10) then 1349 !write(std_out,*) 'Compare bb2 : ' 1350 !write(std_out,*) bb(:), rpt(:), bb(:)-rpt(:) 1351 !write(std_out,*) -cp1232(:)*(n23-sp1232/n23)/n1/n1, -cp3432(:)*sp3432/n23/n2/n2 1352 !end if 1353 bb(:) = rpt(:) 1354 1355 end subroutine ddihedral_d2
ABINIT/deloc2xcart [ Functions ]
NAME
deloc2xcart
FUNCTION
Determine the cartesian coordinates which correspond to the given values of the delocalized coordinates. The relationship is non-linear, so use an iterative scheme, as in Baker JCP .105. 192 (1996). Older reference: Pulay and co. JACS 101 2550 (1979)
INPUTS
deloc <type(delocint)>=Important variables for | pred_delocint | | nang = Number of angles | nbond = Number of bonds | ncart = Number of cartesian directions | (used for constraints) | ndihed = Number of dihedrals | nrshift = Dimension of rshift | ninternal= Number of internal coordinates | ninternal=nbond+nang+ndihed+ncart | | angs(2,3,nang) = Indexes to characterize angles | bonds(2,2,nbond)= For a bond between iatom and jatom | bonds(1,1,nbond) = iatom | bonds(2,1,nbond) = icenter | bonds(1,2,nbond) = jatom | bonds(2,2,nbond) = irshift | carts(2,ncart) = Index of total primitive internal, | and atom (carts(2,:)) | dihedrals(2,4,ndihed)= Indexes to characterize dihedrals | | rshift(3,nrshift)= Shift in xred that must be done to find | all neighbors of a given atom within a | given number of neighboring shells natom = Number of atoms (dtset%natom) rprimd(3,3)=dimensional real space primitive translations (bohr)
OUTPUT
bt_inv_matrix(3*(natom-1),3*natom)=inverse of transpose of B matrix
SIDE EFFECTS
u_matrix(ninternal,3*(natom-1))=eigenvectors of G = BB^T matrix xcart(3,natom)=cartesian coordinates of atoms (bohr)
NOTES
SOURCE
648 subroutine deloc2xcart(deloc,natom,rprimd,xcart,deloc_int,btinv,u_matrix) 649 650 !Arguments ------------------------------------ 651 !scalars 652 integer,intent(in) :: natom 653 type(delocint),intent(in) :: deloc 654 !arrays 655 real(dp),intent(in) :: deloc_int(3*(natom-1)),rprimd(3,3) 656 real(dp),intent(inout) :: u_matrix(deloc%ninternal,3*(natom-1)) 657 real(dp),intent(inout) :: xcart(3,natom) 658 real(dp),intent(out) :: btinv(3*(natom-1),3*natom) 659 660 !Local variables------------------------------- 661 !scalars 662 integer :: iiter,iprim,niter 663 integer :: ii 664 real(dp) :: minmix, maxmix 665 real(dp) :: mix,tot_diff, toldeloc 666 real(dp) :: lntoldeloc 667 logical :: DEBUG=.FALSE. 668 !arrays 669 real(dp) :: btinv_tmp(3*(natom-1),3*natom) 670 real(dp) :: cgrad(3*natom),cgrad_old(3*natom) 671 real(dp) :: deloc_int_now(3*(natom-1)),prim_int(deloc%ninternal) 672 real(dp) :: tmpxcart(3*natom) 673 real(dp) :: xdeloc_diff(3*(natom-1)) 674 675 character(len=500) :: message 676 677 ! ****************************************************************** 678 679 if (DEBUG) then 680 write(ab_out,*) 'ENTERING DELOC2XCART' 681 682 write (message,*) 'BONDS=',deloc%nbond 683 call wrtout(ab_out,message,'COLL') 684 do ii = 1, deloc%nbond 685 write (message,*) ii, deloc%bonds(:,:,ii) 686 call wrtout(ab_out,message,'COLL') 687 end do 688 689 write (message,*) 'ANGS=',deloc%nang 690 call wrtout(ab_out,message,'COLL') 691 do ii = 1, deloc%nang 692 write (message,*) ii, deloc%angs(:,:,ii) 693 call wrtout(ab_out,message,'COLL') 694 end do 695 696 write (message,*) 'DIHEDRALS=',deloc%ndihed 697 call wrtout(ab_out,message,'COLL') 698 do ii = 1, deloc%ndihed 699 write (message,*) ii, deloc%dihedrals(:,:,ii) 700 call wrtout(ab_out,message,'COLL') 701 end do 702 703 write (message,*) 'CARTS=',deloc%ncart 704 call wrtout(ab_out,message,'COLL') 705 do ii = 1, deloc%ncart 706 write (message,*) ii, deloc%carts(:,ii) 707 call wrtout(ab_out,message,'COLL') 708 end do 709 710 write (ab_out,*) 'xcart (input)' 711 do ii=1,natom 712 write (ab_out,*) xcart(:,ii) 713 end do 714 715 end if 716 717 niter = 200 718 tmpxcart = reshape(xcart,(/3*natom/)) 719 720 cgrad_old(:) = zero 721 cgrad(:) = zero 722 maxmix = 0.9_dp 723 minmix = 0.2_dp 724 toldeloc = tol10 725 lntoldeloc = log(toldeloc) 726 727 do iiter=1,niter 728 if (iiter==1) then 729 mix= minmix 730 else 731 mix = minmix + (maxmix-minmix)*(log(tot_diff)-lntoldeloc) / lntoldeloc 732 end if 733 if (mix < minmix) mix = minmix 734 if (mix > maxmix) mix = maxmix 735 736 tmpxcart(:) = tmpxcart(:) + mix*cgrad(:) 737 xcart = reshape(tmpxcart,(/3,natom/)) 738 call xcart2deloc(deloc,natom,rprimd,xcart,& 739 & btinv_tmp,u_matrix,deloc_int_now,prim_int) 740 ! update the BT^{-1} matrix? 741 btinv(:,:) = btinv_tmp(:,:) 742 743 xdeloc_diff(:) = deloc_int(:) - deloc_int_now(:) 744 745 tot_diff = sum(abs(xdeloc_diff)) 746 if (tot_diff < toldeloc) exit 747 748 cgrad_old(:) = cgrad(:) 749 750 ! gradient vector = btinv^{T} * xdeloc_diff 751 call dgemv('T',3*(natom-1),3*natom,one,& 752 & btinv,3*(natom-1),xdeloc_diff,1,zero,cgrad,1) 753 end do 754 !end iiter do 755 756 call xcart2deloc(deloc,natom,rprimd,xcart,& 757 & btinv,u_matrix,deloc_int_now,prim_int) 758 write (message,'(3a)') 'delocalized internals, after convergence of xcart = ', ch10 759 call wrtout(std_out,message,'COLL') 760 do ii = 1, 3*(natom-1) 761 write (message,'(I6,E20.10,2x)') ii, deloc_int_now(ii) 762 call wrtout(std_out,message,'COLL') 763 end do 764 765 xdeloc_diff(:) = deloc_int(:) - deloc_int_now(:) 766 767 write (message,'(a)') 'Primitive internal coordinate values:' 768 call wrtout(std_out,message,'COLL') 769 do iprim = 1, deloc%nbond 770 write (message,'(i6,E20.10)') iprim, prim_int(iprim) 771 call wrtout(std_out,message,'COLL') 772 end do 773 do iprim = deloc%nbond+1, deloc%nbond+deloc%nang+deloc%ndihed 774 write (message,'(i6,2E20.10)') iprim, prim_int(iprim), prim_int(iprim)/pi*180.0_dp 775 call wrtout(std_out,message,'COLL') 776 end do 777 do iprim = deloc%nbond+deloc%nang+deloc%ndihed+1, deloc%ninternal 778 write (message,'(i6,E20.10)') iprim, prim_int(iprim) 779 call wrtout(std_out,message,'COLL') 780 end do 781 782 if (iiter == niter+1) then 783 write (message,'(a,i6,a,E20.10)') 'deloc2xcart : Error, xcart not converged in ', niter, 'iterations ', tot_diff 784 ABI_ERROR(message) 785 end if 786 787 if(DEBUG)then 788 write (ab_out,*) 'xcart (output)' 789 do ii=1,natom 790 write (ab_out,*) xcart(:,ii) 791 end do 792 write(ab_out,*) 'EXITING DELOC2XCART' 793 end if 794 795 end subroutine deloc2xcart
ABINIT/gred2gdeloc [ Functions ]
NAME
gred2gdeloc
FUNCTION
calculate delocalized forces from reduced coordinate ones
INPUTS
btinv(3*(natom-1),3*natom)= inverse transpose of B matrix (see delocint) natom = number of atoms gprimd(3,3)=dimensional translations in reciprocal space (bohr-1)
OUTPUT
deloc_gred(3*(natom-1))=delocalized gradients from reduced coordinate ones gred(3,natom)=delocalized gradients in reduced coordinates
SOURCE
816 subroutine gred2gdeloc(btinv,deloc_gred,gred,natom,gprimd) 817 818 !Arguments ------------------------------------ 819 !scalars 820 integer, intent(in) :: natom 821 !arrays 822 real(dp),intent(in) :: btinv(3*(natom-1),3*natom),gprimd(3,3),gred(3,natom) 823 real(dp),intent(out) :: deloc_gred(3*(natom-1)) 824 825 !Local variables------------------------------- 826 integer :: ii 827 !arrays 828 real(dp) :: fcart(3,natom) 829 character(len=500) :: message 830 831 ! ****************************************************************** 832 833 !make cartesian forces 834 835 call dgemm('N','N',3,natom,3,one,& 836 & gprimd,3,gred,3,zero,fcart,3) 837 838 !turn cartesian to delocalized forces 839 call dgemv('N',3*(natom-1),3*natom,one,& 840 & btinv,3*(natom-1),fcart,1,zero,deloc_gred,1) 841 842 write (message,'(a)') 'gred2gdeloc : deloc_gred = ' 843 call wrtout(std_out,message,'COLL') 844 845 do ii = 1, 3*(natom-1) 846 write (message,'(I6,E16.6)') ii, deloc_gred(ii) 847 call wrtout(std_out,message,'COLL') 848 end do 849 850 end subroutine gred2gdeloc
ABINIT/m_pred_delocint [ Modules ]
NAME
m_pred_delocint
FUNCTION
COPYRIGHT
Copyright (C) 1998-2024 ABINIT group (MVer, DCA, XG, GMR, JCC, SE) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
15 #if defined HAVE_CONFIG_H 16 #include "config.h" 17 #endif 18 19 #include "abi_common.h" 20 21 module m_pred_delocint 22 23 use defs_basis 24 use m_errors 25 use m_abicore 26 use m_abimover 27 use m_abihist 28 use m_xfpack 29 use m_linalg_interfaces 30 31 use m_geometry, only : fcart2gred, xcart2xred, xred2xcart, metric, acrossb 32 use m_bfgs, only : hessinit, hessupdt, brdene 33 use m_results_gs, only : results_gs_type 34 35 implicit none 36 37 private
ABINIT/pred_delocint [ Functions ]
NAME
pred_delocint
FUNCTION
Ionmov predictors (10) BFGS with delocalized internal coordinates IONMOV 10: Given a starting point xred that is a vector of length 3*(natom-1) (reduced nuclei coordinates), and unit cell parameters (acell and rprimd) the Broyden-Fletcher-Goldfarb-Shanno minimization is performed on the total energy function, using its gradient (atomic forces and stresses) while the optimization of unit cell parameters is only performed if optcell/=0. The convergence requirement on the atomic forces, 'tolmxf', allows an early exit. Otherwise no more than 'ntime' steps are performed. Returned quantities are xred, and eventually acell and rprimd (new ones!). Could see Numerical Recipes (Fortran), 1986, page 307. Implements the delocalized internal coordinate scheme of Andzelm et al. in CPL .335. 321 (2001) \ and Baker et al. JCP .105. 192 (1996) B matrix is derivative of delocalized internals wrt cartesian coordinates U matrix is eigenvectors of G = B*B^{T} S matrix is eigenvectors of F = B^{T}B
INPUTS
ab_mover <type(abimover)> : Datatype with all the information needed by the preditor itime : Index of the present iteration ntime : Maximal number of iterations ionmov : (10 or 11) Specific kind of BFGS zDEBUG : if true print some debugging information
OUTPUT
SIDE EFFECTS
hist <type(abihist)> : History of positions,forces acell, rprimd, stresses
SOURCE
89 subroutine pred_delocint(ab_mover,ab_xfh,deloc,forstr,hist,ionmov,itime,zDEBUG,iexit) 90 91 !Arguments ------------------------------------ 92 !scalars 93 type(abimover),intent(in) :: ab_mover 94 type(ab_xfh_type),intent(inout) :: ab_xfh 95 type(abihist),intent(inout) :: hist 96 type(abiforstr),intent(in) :: forstr 97 type(delocint),intent(inout) :: deloc 98 integer,intent(in) :: itime 99 integer,intent(in) :: ionmov 100 integer,intent(in) :: iexit 101 logical,intent(in) :: zDEBUG 102 103 !Local variables------------------------------- 104 !scalars 105 integer :: ndim,cycl_main 106 integer :: ihist_prev,ii,jj,kk 107 real(dp),save :: ucvol0 108 real(dp) :: ucvol 109 real(dp) :: etotal,etotal_prev 110 logical :: DEBUG=.TRUE. 111 !integer,save :: icenter,irshift ! DELOCINT indexes 112 integer,save :: ndeloc ! DELOCINT number of 113 character(len=500) :: message 114 115 !arrays 116 real(dp),allocatable,save :: hessin(:,:),vin(:),vin_prev(:) 117 real(dp),allocatable,save :: vout(:),vout_prev(:) 118 real(dp),save :: acell0(3) ! Initial acell 119 real(dp),save :: rprimd0(3,3) ! Initial rprimd 120 real(dp),allocatable :: prim_int(:) 121 real(dp),allocatable,save :: u_matrix(:,:) ! DELOCINT this may need to be added to type inside ab_mover 122 real(dp) :: acell(3) 123 real(dp) :: rprimd(3,3) 124 real(dp) :: gprimd(3,3) 125 real(dp) :: gmet(3,3) 126 real(dp) :: rmet(3,3) 127 real(dp) :: residual(3,ab_mover%natom) 128 !real(dp) :: residual_corrected(3,ab_mover%natom) 129 real(dp) :: xred(3,ab_mover%natom),xcart(3,ab_mover%natom) 130 real(dp) :: strten(6) 131 real(dp) :: deloc_gred(3*(ab_mover%natom-1)) 132 real(dp) :: deloc_int(3*(ab_mover%natom-1)) 133 real(dp) :: bt_inv_matrix(3*(ab_mover%natom-1),3*ab_mover%natom) 134 135 !*************************************************************************** 136 !Beginning of executable session 137 !*************************************************************************** 138 139 if(iexit/=0)then 140 ABI_SFREE(vin) 141 ABI_SFREE(vout) 142 ABI_SFREE(vin_prev) 143 ABI_SFREE(vout_prev) 144 ABI_SFREE(hessin) 145 ABI_SFREE(u_matrix) 146 return 147 end if 148 149 !write(std_out,*) 'delocint 01' 150 !########################################################## 151 !### 01. Debugging and Verbose 152 153 if(DEBUG)then 154 write(std_out,'(a,3a,38a,39a)') ch10,('-',kk=1,3),'Debugging and Verbose for pred_deloint',('-',kk=1,39) 155 write(std_out,*) 'ionmov: ',ionmov 156 write(std_out,*) 'itime: ',itime 157 end if 158 159 !write(std_out,*) 'delocint 02' 160 !########################################################## 161 !### 02. Compute the dimension of vectors (ndim) 162 163 !With internal we have 1 coordinate less 164 ndeloc = 3*(ab_mover%natom-1) 165 ndim=ndeloc 166 deloc_int(:)=zero 167 deloc_gred(:)=zero 168 if(ab_mover%optcell==1) ndim=ndim+1 169 if(ab_mover%optcell==2 .or.& 170 & ab_mover%optcell==3) ndim=ndim+6 171 if(ab_mover%optcell>=4) ndim=ndim+3 172 173 if(DEBUG) write(std_out,*) 'Dimension of vin, vout and hessian (ndim): ',ndim 174 175 !write(std_out,*) 'delocint 03' 176 !########################################################## 177 !### 03. Allocate the vectors vin, vout and hessian matrix 178 179 !Notice thqt vin, vout, etc could be allocated 180 !From a previous dataset with a different ndim 181 if(itime==1)then 182 ABI_SFREE(vin) 183 ABI_SFREE(vout) 184 ABI_SFREE(vin_prev) 185 ABI_SFREE(vout_prev) 186 ABI_SFREE(hessin) 187 188 ABI_MALLOC(vin,(ndim)) 189 ABI_MALLOC(vout,(ndim)) 190 ABI_MALLOC(vin_prev,(ndim)) 191 ABI_MALLOC(vout_prev,(ndim)) 192 ABI_MALLOC(hessin,(ndim,ndim)) 193 end if 194 195 196 !write(std_out,*) 'delocint 04' 197 !########################################################## 198 !### 04. Obtain the present values from the history 199 200 call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 201 call xred2xcart(ab_mover%natom,rprimd,xcart,xred) 202 203 strten(:)=hist%strten(:,hist%ihist) 204 etotal =hist%etot(hist%ihist) 205 206 !Fill the residual with forces (No preconditioning) 207 !Or the preconditioned forces 208 if (ab_mover%goprecon==0)then 209 call fcart2gred(hist%fcart(:,:,hist%ihist),residual,rprimd,ab_mover%natom) 210 else 211 residual(:,:)= forstr%gred(:,:) 212 end if 213 214 if(zDEBUG)then 215 write (std_out,*) 'residual:' 216 do kk=1,ab_mover%natom 217 write (std_out,*) residual(:,kk) 218 end do 219 write (std_out,*) 'strten:' 220 write (std_out,*) strten(1:3),ch10,strten(4:6) 221 write (std_out,*) 'etotal:' 222 write (std_out,*) etotal 223 end if 224 225 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 226 227 !Save initial values 228 if (itime==1)then 229 acell0(:)=acell(:) 230 rprimd0(:,:)=rprimd(:,:) 231 ucvol0=ucvol 232 end if 233 234 !DEBUG (UCVOL) 235 if(DEBUG)then 236 write(std_out,*) 'Volume of cell (ucvol):',ucvol 237 end if 238 239 !Get rid of mean force on whole unit cell, but only if no 240 !generalized constraints are in effect 241 ! residual_corrected(:,:)=residual(:,:) 242 ! if(ab_mover%nconeq==0)then 243 ! do ii=1,3 244 ! if (ii/=3.or.ab_mover%jellslab==0) then 245 ! favg=sum(residual_corrected(ii,:))/dble(ab_mover%natom) 246 ! residual_corrected(ii,:)=residual_corrected(ii,:)-favg 247 ! end if 248 ! end do 249 ! end if 250 251 !write(std_out,*) 'delocint 05' 252 !########################################################## 253 !### 05. Compute internals for first time 254 255 if (itime==1)then 256 call make_prim_internals(deloc,ab_mover%natom,& 257 & ab_mover%ntypat,rprimd,ab_mover%typat,xcart,ab_mover%znucl) 258 259 ABI_MALLOC(prim_int,(deloc%ninternal)) 260 261 if(DEBUG)then 262 write (message,'(a,i6)') 'Number of primitive internal coordinates (ninternal): ',deloc%ninternal 263 call wrtout(std_out, message,'COLL') 264 end if 265 266 if (allocated(u_matrix)) then 267 ABI_FREE(u_matrix) 268 end if 269 ABI_MALLOC(u_matrix,(deloc%ninternal,ndeloc)) 270 271 call calc_prim_int(deloc,ab_mover%natom,rprimd,xcart,prim_int) 272 273 if(DEBUG)then 274 write (message,'(a)') 'Primitive internal coordinate values:' 275 call wrtout(std_out, message,'COLL') 276 write (message,'(a)') ' Bonds:' 277 call wrtout(std_out, message,'COLL') 278 do ii = 1, deloc%nbond 279 write (message,'(i6,E20.10)') ii, prim_int(ii) 280 call wrtout(std_out, message,'COLL') 281 end do 282 283 write (message,'(a)') ' Angles:' 284 call wrtout(std_out, message,'COLL') 285 do ii = deloc%nbond+1, deloc%nbond+deloc%nang 286 write (message,'(i6,2(E20.10,2x))') ii, prim_int(ii), prim_int(ii)/pi*180.0_dp 287 call wrtout(std_out, message,'COLL') 288 end do 289 290 write (message,'(a)') ' Dihedrals:' 291 call wrtout(std_out, message,'COLL') 292 do ii = deloc%nbond+deloc%nang+1, deloc%nbond+deloc%nang+deloc%ndihed 293 write (message,'(i6,2(E20.10,2x))') ii, prim_int(ii), prim_int(ii)/pi*180.0_dp 294 call wrtout(std_out, message,'COLL') 295 end do 296 297 write (message,'(a)') ' Cartesian auxiliary coordinates for constraints:' 298 call wrtout(std_out, message,'COLL') 299 do ii = deloc%nbond+deloc%nang+deloc%ndihed+1, deloc%ninternal 300 write (message,'(i6,E20.10)') ii, prim_int(ii) 301 call wrtout(std_out, message,'COLL') 302 end do 303 end if 304 305 ABI_FREE(prim_int) 306 307 ! equal weight on all internal coordinates as a starting point. 308 u_matrix(:,:) = one / dble (ndeloc) 309 310 ! Zero the arrays before first use 311 deloc_gred(:) = zero 312 313 end if 314 315 ABI_MALLOC(prim_int,(deloc%ninternal)) 316 317 !write(std_out,*) 'delocint 06' 318 !########################################################## 319 !### 06. Compute delocalized coordinates and forces 320 321 !xcart ---> deloc_int 322 323 !Convert positions to delocalized coordinates for next step 324 call xcart2deloc(deloc,ab_mover%natom,rprimd,xcart,& 325 & bt_inv_matrix,u_matrix,deloc_int,prim_int) 326 327 !gred ---> deloc_gred 328 329 !Convert gradients to delocalized coordinates for next step 330 call gred2gdeloc(bt_inv_matrix,deloc_gred,residual,ab_mover%natom,gprimd) 331 332 !write(std_out,*) 'delocint 07' 333 !########################################################## 334 !### 07. Fill the vectors vin and vout 335 336 !DEBUG deloc_int and deloc_gred before pack 337 if(DEBUG)then 338 write (std_out,*) 'Delocalized internals and forces (ndeloc):',ndeloc 339 write(std_out,*) 'deloc_int' 340 do ii=1,ndeloc,3 341 if (ii+2<=ndeloc)then 342 write(std_out,*) ii,deloc_int(ii:ii+2) 343 else 344 write(std_out,*) ii,deloc_int(ii:ndeloc) 345 end if 346 end do 347 write(std_out,*) 'deloc_gred' 348 do ii=1,ndeloc,3 349 if (ii+2<=ndeloc)then 350 write(std_out,*) ii,deloc_gred(ii:ii+2) 351 else 352 write(std_out,*) ii,deloc_gred(ii:ndeloc) 353 end if 354 end do 355 end if 356 357 !DELOCINT 358 !Instead of gred_corrected we use deloc_gred 359 !Instead of xred e use deloc_int 360 ! 361 !Initialize input vectors : first vin, then vout 362 !The values of vin from the previous iteration 363 !should be the same 364 call xfpack_x2vin(acell, ab_mover%natom-1, ndim,& 365 & ab_mover%nsym, ab_mover%optcell, rprimd, rprimd0,& 366 & ab_mover%symrel, ucvol, ucvol0, vin, deloc_int) 367 !end if 368 369 call xfpack_f2vout(deloc_gred, ab_mover%natom-1, ndim,& 370 & ab_mover%optcell, ab_mover%strtarget, strten, ucvol, & 371 & vout) 372 373 !write(std_out,*) 'delocint 08' 374 !########################################################## 375 !### 08. Initialize or update the hessian matrix 376 377 !Initialise the Hessian matrix using gmet 378 if (itime==1)then 379 380 ! Initialise the Hessian matrix with ab_mover%userrc. 381 ! this has become unusable because it imposes ndim >= 3 natom 382 ! ident = 3x3 identity matrix 383 ! call hessinit(ab_mover, hessin, gmet, ndim, ucvol) 384 hessin = zero 385 do ii=1, ndim 386 hessin (ii,ii) = one 387 end do 388 389 ! ! Initialize inverse hessian with identity matrix 390 ! ! in cartesian coordinates, which makes use of metric tensor gmet 391 ! ! in reduced coordinates. 392 ! hessin(:,:)=zero 393 ! do ii=1,ab_mover%natom 394 ! do kk=1,3 395 ! do jj=1,3 396 ! ! Warning : implemented in reduced coordinates 397 ! if (ab_mover%iatfix(kk,ii)==0 .and.& 398 ! & ab_mover%iatfix(jj,ii)==0 )then 399 ! hessin(kk+3*(ii-1),jj+3*(ii-1))=gmet(kk,jj) 400 ! end if 401 ! end do 402 ! end do 403 ! end do 404 ! if(ab_mover%optcell/=0)then 405 ! ! These values might lead to too large changes in some cases 406 ! diag=ab_mover%strprecon*30.0_dp/ucvol 407 ! if(ab_mover%optcell==1) diag=diag/three 408 ! do ii=3*ab_mover%natom+1,ndim 409 ! hessin(ii,ii)=diag 410 ! end do 411 ! end if 412 413 if (ab_mover%restartxf/=0) then 414 415 call xfh_recover_deloc(ab_xfh,ab_mover,acell,cycl_main,& 416 & residual,hessin,ndim,rprimd,rprimd0,strten,ucvol,ucvol0,vin,vin_prev,& 417 & vout,vout_prev,xred,deloc,deloc_int,deloc_gred,bt_inv_matrix,gprimd,prim_int,& 418 & u_matrix) 419 420 end if 421 422 end if 423 424 ABI_FREE(prim_int) 425 426 if(itime>1)then 427 ! Update the hessian matrix, by taking into account the 428 ! current pair (x,f) and the previous one. 429 call hessupdt(hessin,ab_mover%iatfix,ab_mover%natom,ndim,vin,& 430 & vin_prev,vout,vout_prev) 431 432 end if 433 434 !DEBUG (vin,vout and hessin before prediction) 435 if(DEBUG)then 436 write(std_out,*) 'Vectors vin and vout and inverse of Hessian (hessin) [before]' 437 write(std_out,*) 'vin:' 438 do ii=1,ndim,3 439 if (ii+2<=ndim)then 440 write(std_out,*) ii,vin(ii:ii+2) 441 else 442 write(std_out,*) ii,vin(ii:ndim) 443 end if 444 end do 445 write(std_out,*) 'vout:' 446 do ii=1,ndim,3 447 if (ii+2<=ndim)then 448 write(std_out,*) ii,vout(ii:ii+2) 449 else 450 write(std_out,*) ii,vout(ii:ndim) 451 end if 452 end do 453 write(std_out,*) 'Inverse Hessian (hessin): ',ndim,'x',ndim 454 do kk=1,ndim 455 do jj=1,ndim,3 456 if (jj+2<=ndim)then 457 write(std_out,*) jj,hessin(jj:jj+2,kk) 458 else 459 write(std_out,*) jj,hessin(jj:ndim,kk) 460 end if 461 end do 462 end do 463 end if 464 465 !write(std_out,*) 'delocint 09' 466 !########################################################## 467 !### 09. Compute the next values 468 469 if(ionmov==10 .or. itime==1)then 470 471 ! Previous cartesian coordinates 472 vin_prev(:)=vin(:) 473 474 ! New atomic cartesian coordinates are obtained from vin, hessin 475 ! and vout 476 do ii=1,ndim 477 vin(:)=vin(:)-hessin(:,ii)*vout(ii) 478 end do 479 ! Previous atomic forces 480 vout_prev(:)=vout(:) 481 482 else 483 if(ionmov==11)then 484 ihist_prev = abihist_findIndex(hist,-1) 485 etotal_prev=hist%etot(ihist_prev) 486 ! Here the BFGS algorithm, modified to take into account the 487 ! energy 488 call brdene(etotal,etotal_prev,hessin,& 489 & ndim,vin,vin_prev,vout,vout_prev) 490 491 end if 492 493 ! DEBUG (vin,vout and hessin after prediction) 494 if(DEBUG)then 495 write(std_out,*) 'Vectors vin and vout [after prediction]' 496 write(std_out,*) 'vin:' 497 do ii=1,ndim,3 498 if (ii+2<=ndim)then 499 write(std_out,*) ii,vin(ii:ii+2) 500 else 501 write(std_out,*) ii,vin(ii:ndim) 502 end if 503 end do 504 write(std_out,*) 'vout:' 505 do ii=1,ndim,3 506 if (ii+2<=ndim)then 507 write(std_out,*) ii,vout(ii:ii+2) 508 else 509 write(std_out,*) ii,vout(ii:ndim) 510 end if 511 end do 512 end if 513 514 ! Implement fixing of atoms : put back old values for fixed 515 ! components 516 do kk=1,ab_mover%natom 517 do jj=1,3 518 ! Warning : implemented in reduced coordinates 519 if ( ab_mover%iatfix(jj,kk)==1) then 520 vin(jj+(kk-1)*3)=vin_prev(jj+(kk-1)*3) 521 end if 522 end do 523 end do 524 end if 525 526 !write(std_out,*) 'delocint 10' 527 !########################################################## 528 !### 10. Convert from delocalized to xcart and xred 529 530 !Transfer vin to deloc_int, acell and rprimd 531 call xfpack_vin2x(acell, acell0, ab_mover%natom-1, ndim,& 532 & ab_mover%nsym, ab_mover%optcell, rprimd, rprimd0,& 533 & ab_mover%symrel, ucvol, ucvol0,& 534 & vin, deloc_int) 535 536 if(DEBUG)then 537 write (std_out,*) 'Delocalized internals (deloc_int) [after prediction]:' 538 write(std_out,*) 'deloc_int:' 539 do ii=1,ndeloc,3 540 if (ii+2<=ndeloc)then 541 write(std_out,*) ii,deloc_int(ii:ii+2) 542 else 543 write(std_out,*) ii,deloc_int(ii:ndeloc) 544 end if 545 end do 546 write(std_out,*) 'BT Inverse Matrix:' 547 do ii=1,3*ab_mover%natom 548 write(std_out,*) bt_inv_matrix(:,ii) 549 end do 550 write (std_out,*) 'xcart (before deloc2xcart):' 551 do ii=1,ab_mover%natom 552 write (std_out,*) xcart(:,ii) 553 end do 554 end if 555 556 !this routine contains an iterative scheme to find xcart 557 !from the non-linear relations between deloc and xcart 558 !SIGNIFICANTLY DIFFERENT FROM xcart2deloc 559 call deloc2xcart(deloc,ab_mover%natom,rprimd,xcart,& 560 & deloc_int,bt_inv_matrix,u_matrix) 561 562 !Convert new xcart (cartesian) to xred (reduced coordinates) 563 call xcart2xred(ab_mover%natom,rprimd,xcart,xred) 564 565 566 !write(std_out,*) 'delocint 11' 567 !########################################################## 568 !### 11. Update the history with the prediction 569 570 !Increase indexes 571 hist%ihist = abihist_findIndex(hist,+1) 572 573 if(ab_mover%optcell/=0)then 574 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 575 end if 576 577 !Fill the history with the variables 578 !xred, acell, rprimd, vel 579 call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 580 ihist_prev = abihist_findIndex(hist,-1) 581 hist%vel(:,:,hist%ihist)=hist%vel(:,:,ihist_prev) 582 583 if(zDEBUG)then 584 write (std_out,*) 'residual:' 585 do kk=1,ab_mover%natom 586 write (std_out,*) residual(:,kk) 587 end do 588 write (std_out,*) 'strten:' 589 write (std_out,*) strten(1:3),ch10,strten(4:6) 590 write (std_out,*) 'etotal:' 591 write (std_out,*) etotal 592 end if 593 594 end subroutine pred_delocint
ABINIT/xcart2deloc [ Functions ]
NAME
xcart2deloc
FUNCTION
Calculate values of delocalized coordinates as a function of cartesian ones. First primitive internals, then B matrix, then F, then U then delocalized internals.
INPUTS
deloc <type(delocint)>=Important variables for pred_delocint | | nang = Number of angles | nbond = Number of bonds | ncart = Number of cartesian directions | (used for constraints) | ndihed = Number of dihedrals | nrshift = Dimension of rshift | ninternal= Number of internal coordinates | ninternal=nbond+nang+ndihed+ncart | | angs(2,3,nang) = Indexes to characterize angles | bonds(2,2,nbond)= For a bond between iatom and jatom | bonds(1,1,nbond) = iatom | bonds(2,1,nbond) = icenter | bonds(1,2,nbond) = jatom | bonds(2,2,nbond) = irshift | carts(2,ncart) = Index of total primitive internal, | and atom (carts(2,:)) | dihedrals(2,4,ndihed)= Indexes to characterize dihedrals | | rshift(3,nrshift)= Shift in xred that must be done to find | all neighbors of a given atom within a | given number of neighboring shells natom = Number of atoms rprimd(3,3) = Dimensional real space primitive translations (bohr) xcart(3,natom) = Cartesian coordinates of atoms (bohr)
OUTPUT
bt_inv_matrix(3*(natom-1),3*natom) = Inverse of B^{T} matrix
deloc_int(3*(natom-1)) = Delocalized internal coordinates
prim_int(ninternal) = Primitive internal coordinates
SIDE EFFECTS
u_matrix(ninternal,3*(natom-1)) = Eigenvectors of BB^T matrix
NOTES
SOURCE
1409 subroutine xcart2deloc(deloc,natom,rprimd,xcart,bt_inv_matrix,u_matrix,deloc_int,prim_int) 1410 1411 !Arguments ------------------------------------ 1412 !scalars 1413 integer,intent(in) :: natom 1414 type(delocint),intent(in) :: deloc 1415 !arrays 1416 real(dp),intent(in) :: rprimd(3,3),xcart(3,natom) 1417 real(dp),intent(inout) :: u_matrix(deloc%ninternal,3*(natom-1)) 1418 real(dp),intent(out) :: bt_inv_matrix(3*(natom-1),3*natom) 1419 real(dp),intent(out) :: deloc_int(3*(natom-1)) 1420 real(dp),intent(out) :: prim_int(deloc%ninternal) 1421 1422 !Local variables------------------------------- 1423 !scalars 1424 integer :: ii 1425 logical :: DEBUG=.FALSE. 1426 !arrays 1427 real(dp) :: b_matrix(deloc%ninternal,3*natom) 1428 1429 ! ****************************************************************** 1430 1431 call calc_prim_int(deloc,natom,rprimd,xcart,prim_int) 1432 if (DEBUG)then 1433 write(std_out,*) 'Primitive Internals' 1434 do ii=1,deloc%ninternal 1435 write(std_out,*) prim_int(ii) 1436 end do 1437 end if 1438 1439 call calc_b_matrix(deloc,natom,rprimd,xcart,b_matrix) 1440 if (DEBUG)then 1441 write(std_out,*) 'B Matrix' 1442 do ii=1,deloc%ninternal 1443 write(std_out,*) b_matrix(:,ii) 1444 end do 1445 end if 1446 1447 call calc_btinv_matrix(b_matrix,natom,deloc%ninternal,& 1448 & bt_inv_matrix,u_matrix) 1449 if (DEBUG)then 1450 write(std_out,*) 'BT Inverse Matrix' 1451 do ii=1,3*natom 1452 write(std_out,*) bt_inv_matrix(:,ii) 1453 end do 1454 end if 1455 1456 !calculate value of delocalized internals 1457 1458 call dgemv('T',deloc%ninternal,3*(natom-1),one,& 1459 & u_matrix,deloc%ninternal,prim_int,1,zero,deloc_int,1) 1460 1461 end subroutine xcart2deloc
ABINIT/xfh_recover_deloc [ Functions ]
NAME
xfh_recover_deloc
FUNCTION
Update the contents of the history xfhist taking values from xred, acell, rprim, gred_corrected and strten
INPUTS
OUTPUT
SOURCE
1649 subroutine xfh_recover_deloc(ab_xfh,ab_mover,acell,cycl_main,& 1650 & gred,hessin,ndim,rprim,rprimd0,strten,ucvol,ucvol0,vin,vin_prev,& 1651 & vout,vout_prev,xred,deloc,deloc_int,deloc_gred,btinv,gprimd,prim_int,& 1652 & u_matrix) 1653 1654 !Arguments ------------------------------------ 1655 !scalars 1656 1657 integer,intent(in) :: ndim 1658 integer,intent(out) :: cycl_main 1659 real(dp),intent(inout) :: ucvol,ucvol0 1660 type(ab_xfh_type),intent(inout) :: ab_xfh 1661 type(abimover),intent(in) :: ab_mover 1662 ! DELOCINT specials 1663 type(delocint),intent(in) :: deloc 1664 1665 !arrays 1666 real(dp),intent(inout) :: acell(3) 1667 real(dp),intent(inout) :: hessin(:,:) 1668 real(dp),intent(inout) :: xred(3,ab_mover%natom) 1669 real(dp),intent(inout) :: rprim(3,3) 1670 real(dp),intent(inout) :: rprimd0(3,3) 1671 real(dp),intent(inout) :: gred(3,ab_mover%natom) 1672 real(dp),intent(inout) :: strten(6) 1673 real(dp),intent(inout) :: vin(:) 1674 real(dp),intent(inout) :: vin_prev(:) 1675 real(dp),intent(inout) :: vout(:) 1676 real(dp),intent(inout) :: vout_prev(:) 1677 ! DELOCINT specials 1678 real(dp),intent(inout) :: deloc_gred(3*(ab_mover%natom-1)) 1679 real(dp),intent(inout) :: deloc_int(3*(ab_mover%natom-1)) 1680 real(dp),intent(inout) :: btinv(3*(ab_mover%natom-1),3*ab_mover%natom) 1681 real(dp),intent(inout) :: prim_int(:),u_matrix(:,:),gprimd(3,3) 1682 1683 !Local variables------------------------------- 1684 !scalars 1685 integer :: ixfh 1686 real(dp) :: xcart(3,ab_mover%natom) 1687 1688 !********************************************************************* 1689 1690 if(ab_xfh%nxfh/=0)then 1691 ! Loop over previous time steps 1692 do ixfh=1,ab_xfh%nxfh 1693 1694 ! For that time step, get new (x,f) from xfhist 1695 xred(:,:) =ab_xfh%xfhist(:,1:ab_mover%natom ,1,ixfh) 1696 rprim(1:3,1:3)=ab_xfh%xfhist(:,ab_mover%natom+2:ab_mover%natom+4,1,ixfh) 1697 acell(:) =ab_xfh%xfhist(:,ab_mover%natom+1,1,ixfh) 1698 gred(:,:) =ab_xfh%xfhist(:,1:ab_mover%natom,2,ixfh) 1699 ! This use of results_gs is unusual 1700 strten(1:3) =ab_xfh%xfhist(:,ab_mover%natom+2,2,ixfh) 1701 strten(4:6) =ab_xfh%xfhist(:,ab_mover%natom+3,2,ixfh) 1702 1703 ! !DEBUG 1704 ! write (ab_out,*) '---READ FROM XFHIST---' 1705 1706 ! write (ab_out,*) 'XRED' 1707 ! do kk=1,ab_mover%natom 1708 ! write (ab_out,*) xred(:,kk) 1709 ! end do 1710 ! write (ab_out,*) 'FRED' 1711 ! do kk=1,ab_mover%natom 1712 ! write (ab_out,*) gred(:,kk) 1713 ! end do 1714 ! write(ab_out,*) 'RPRIM' 1715 ! do kk=1,3 1716 ! write(ab_out,*) rprim(:,kk) 1717 ! end do 1718 ! write(ab_out,*) 'ACELL' 1719 ! write(ab_out,*) acell(:) 1720 ! !DEBUG 1721 1722 ! Convert input xred (reduced coordinates) to xcart (cartesian) 1723 call xred2xcart(ab_mover%natom,rprimd0,xcart,xred) 1724 ! Convert input coordinates in Delocalized internals 1725 call xcart2deloc(deloc,ab_mover%natom,rprimd0,xcart,& 1726 & btinv,u_matrix,deloc_int,prim_int) 1727 ! Convert forces to delocalized coordinates for next step 1728 call gred2gdeloc(btinv,deloc_gred,gred,ab_mover%natom,gprimd) 1729 1730 ! Transfer it in vin, vout 1731 call xfpack_x2vin(acell,ab_mover%natom-1,& 1732 & ndim,ab_mover%nsym,ab_mover%optcell,rprim,rprimd0,& 1733 & ab_mover%symrel,ucvol,ucvol0,vin,deloc_int) 1734 call xfpack_f2vout(deloc_gred,ab_mover%natom-1,& 1735 & ndim,ab_mover%optcell,ab_mover%strtarget,strten,& 1736 & ucvol,vout) 1737 ! Get old time step, if any, and update inverse hessian 1738 if(ixfh/=1)then 1739 xred(:,:) =ab_xfh%xfhist(:,1:ab_mover%natom,1,ixfh-1) 1740 rprim(1:3,1:3)=& 1741 & ab_xfh%xfhist(:,ab_mover%natom+2:ab_mover%natom+4,1,ixfh-1) 1742 acell(:)=ab_xfh%xfhist(:,ab_mover%natom+1,1,ixfh-1) 1743 gred(:,:)=ab_xfh%xfhist(:,1:ab_mover%natom,2,ixfh-1) 1744 ! This use of results_gs is unusual 1745 strten(1:3)=ab_xfh%xfhist(:,ab_mover%natom+2,2,ixfh-1) 1746 strten(4:6)=ab_xfh%xfhist(:,ab_mover%natom+3,2,ixfh-1) 1747 1748 ! Convert input xred (reduced coordinates) to xcart (cartesian) 1749 call xred2xcart(ab_mover%natom,rprimd0,xcart,xred) 1750 ! Convert input coordinates in Delocalized internals 1751 call xcart2deloc(deloc,ab_mover%natom,rprimd0,xcart,& 1752 & btinv,u_matrix,deloc_int,prim_int) 1753 ! Convert forces to delocalized coordinates for next step 1754 call gred2gdeloc(btinv,deloc_gred,gred,ab_mover%natom,gprimd) 1755 1756 ! Tranfer it in vin_prev, vout_prev 1757 call xfpack_x2vin(acell,ab_mover%natom-1,& 1758 & ndim,ab_mover%nsym,ab_mover%optcell,rprim,rprimd0,& 1759 & ab_mover%symrel,ucvol,ucvol0,vin_prev,deloc_int) 1760 call xfpack_f2vout(deloc_gred,ab_mover%natom-1,& 1761 & ndim,ab_mover%optcell,ab_mover%strtarget,strten,& 1762 & ucvol,vout_prev) 1763 1764 ! write(ab_out,*) 'Hessian matrix before update',ndim,'x',ndim 1765 ! write(ab_out,*) 'ixfh=',ixfh 1766 ! do kk=1,ndim 1767 ! do jj=1,ndim,3 1768 ! if (jj+2<=ndim)then 1769 ! write(ab_out,*) jj,hessin(jj:jj+2,kk) 1770 ! else 1771 ! write(ab_out,*) jj,hessin(jj:ndim,kk) 1772 ! end if 1773 ! end do 1774 ! end do 1775 1776 call hessupdt(hessin,ab_mover%iatfix,ab_mover%natom-1,ndim,& 1777 & vin,vin_prev,vout,vout_prev) 1778 1779 ! !DEBUG 1780 ! write(ab_out,*) 'Hessian matrix after update',ndim,'x',ndim 1781 ! do kk=1,ndim 1782 ! do jj=1,ndim,3 1783 ! if (jj+2<=ndim)then 1784 ! write(ab_out,*) jj,hessin(jj:jj+2,kk) 1785 ! else 1786 ! write(ab_out,*) jj,hessin(jj:ndim,kk) 1787 ! end if 1788 ! end do 1789 ! end do 1790 ! !DEBUG 1791 1792 end if !if(ab_xfh%nxfh/=0) 1793 1794 ! End loop over previous time steps 1795 end do 1796 1797 ! The hessian has been generated, 1798 ! as well as the latest vin and vout 1799 ! so will cycle the main loop 1800 cycl_main=1 1801 1802 end if 1803 1804 end subroutine xfh_recover_deloc