TABLE OF CONTENTS


ABINIT/align_u_matrices [ Functions ]

[ Top ] [ Functions ]

NAME

 align_u_matrices

FUNCTION

INPUTS

OUTPUT

PARENTS

      xcart2deloc

CHILDREN

SOURCE

1774  subroutine align_u_matrices(natom,ninternal,u_matrix,u_matrix_old,s_matrix,f_eigs)
1775 
1776 
1777 !This section has been created automatically by the script Abilint (TD).
1778 !Do not modify the following lines by hand.
1779 #undef ABI_FUNC
1780 #define ABI_FUNC 'align_u_matrices'
1781 !End of the abilint section
1782 
1783  implicit none
1784 
1785 !Arguments ------------------------------------
1786 !scalars
1787  integer,intent(in) :: ninternal,natom
1788 !arrays
1789  real(dp),intent(in) :: u_matrix_old(ninternal,3*(natom-1))
1790  real(dp),intent(inout) :: f_eigs(3*natom)
1791  real(dp),intent(inout) :: s_matrix(3*natom,3*natom)
1792  real(dp),intent(inout) :: u_matrix(ninternal,3*(natom-1))
1793 
1794 !Local variables ------------------------------
1795 !scalars
1796  integer :: ii,iint1,imax
1797  real(dp) :: ss
1798 !arrays
1799  integer :: eigv_flag(3*(natom-1)),eigv_ind(3*(natom-1))
1800  real(dp) :: tmps(3*natom,3*natom)
1801  real(dp) :: tmpu(ninternal,3*(natom-1))
1802  real(dp) :: tmpf(3*natom)
1803 
1804 !******************************************************************
1805 
1806  eigv_flag(:) = 0
1807  eigv_ind(:) = 0
1808 
1809 !just permit a change in sign
1810  do iint1=1,3*(natom-1)
1811    ss = zero
1812    do ii=1,ninternal
1813      ss = ss + u_matrix_old(ii,iint1)*u_matrix(ii,iint1)
1814    end do
1815    if (ss < -tol12) then
1816      imax = -iint1
1817    else
1818      imax = iint1
1819    end if
1820    eigv_ind(iint1) = imax
1821    eigv_flag(abs(imax)) = 1
1822  end do
1823 
1824  tmpu(:,:) = u_matrix
1825  tmps(:,:) = s_matrix
1826  tmpf(:) = f_eigs
1827 !exchange eigenvectors...
1828  do iint1=1,3*(natom-1)
1829    ss = one
1830    if (eigv_ind(iint1) < 0) ss = -one
1831 
1832    imax = abs(eigv_ind(iint1))
1833 
1834    tmpu(:,imax) = ss*u_matrix(:,iint1)
1835 
1836    tmps(:,imax+3) = ss*s_matrix(:,iint1+3)
1837 
1838    tmpf(imax+3) = f_eigs(iint1+3)
1839  end do
1840 
1841  u_matrix(:,:) = tmpu(:,:)
1842  s_matrix(:,:) = tmps(:,:)
1843  f_eigs(:) = tmpf(:)
1844 
1845 end subroutine align_u_matrices

ABINIT/calc_b_matrix [ Functions ]

[ Top ] [ 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

PARENTS

      xcart2deloc

CHILDREN

      acrossb

SOURCE

 974 subroutine calc_b_matrix(deloc,natom,rprimd,xcart,b_matrix)
 975 
 976 
 977 !This section has been created automatically by the script Abilint (TD).
 978 !Do not modify the following lines by hand.
 979 #undef ABI_FUNC
 980 #define ABI_FUNC 'calc_b_matrix'
 981 !End of the abilint section
 982 
 983  implicit none
 984 
 985 !Arguments ------------------------------------
 986 !scalars
 987  integer,intent(in) :: natom
 988  type(delocint),intent(in) :: deloc
 989 
 990 !arrays
 991  real(dp),intent(in) :: rprimd(3,3),xcart(3,natom)
 992  real(dp),intent(out) :: b_matrix(deloc%ninternal,3*natom)
 993 
 994 !Local variables-------------------------------
 995 !scalars
 996  integer :: i1,i2,i3,i4,iang,ibond,icart,idihed,iprim,s1,s2,s3,s4
 997 !arrays
 998  real(dp) :: bb(3),r1(3),r2(3),r3(3),r4(3)
 999 
1000 ! *************************************************************************
1001 
1002  iprim=0
1003  b_matrix(:,:) = zero
1004 
1005  do ibond=1,deloc%nbond
1006    i1 = deloc%bonds(1,1,ibond)
1007    s1 = deloc%bonds(2,1,ibond)
1008    r1(:) = xcart(:,i1)+deloc%rshift(1,s1)*rprimd(:,1)&
1009 &   +deloc%rshift(2,s1)*rprimd(:,2)&
1010 &   +deloc%rshift(3,s1)*rprimd(:,3)
1011    i2 = deloc%bonds(1,2,ibond)
1012    s2 = deloc%bonds(2,2,ibond)
1013    r2(:) = xcart(:,i2)+deloc%rshift(1,s2)*rprimd(:,1)&
1014 &   +deloc%rshift(2,s2)*rprimd(:,2)&
1015 &   +deloc%rshift(3,s2)*rprimd(:,3)
1016    iprim=iprim+1
1017    call dbond_length_d1(r1,r2,bb)
1018    b_matrix(iprim,3*(i1-1)+1:3*i1) = b_matrix(iprim,3*(i1-1)+1:3*i1) + bb(:)
1019    call dbond_length_d1(r2,r1,bb)
1020    b_matrix(iprim,3*(i2-1)+1:3*i2) = b_matrix(iprim,3*(i2-1)+1:3*i2) + bb(:)
1021  end do
1022 
1023 !second: angle values (ang)
1024  do iang=1,deloc%nang
1025    i1 = deloc%angs(1,1,iang)
1026    s1 = deloc%angs(2,1,iang)
1027    r1(:) = xcart(:,i1)+deloc%rshift(1,s1)*rprimd(:,1)&
1028 &   +deloc%rshift(2,s1)*rprimd(:,2)&
1029 &   +deloc%rshift(3,s1)*rprimd(:,3)
1030    i2 = deloc%angs(1,2,iang)
1031    s2 = deloc%angs(2,2,iang)
1032    r2(:) = xcart(:,i2)+deloc%rshift(1,s2)*rprimd(:,1)&
1033 &   +deloc%rshift(2,s2)*rprimd(:,2)&
1034 &   +deloc%rshift(3,s2)*rprimd(:,3)
1035    i3 = deloc%angs(1,3,iang)
1036    s3 = deloc%angs(2,3,iang)
1037    r3(:) = xcart(:,i3)+deloc%rshift(1,s3)*rprimd(:,1)&
1038 &   +deloc%rshift(2,s3)*rprimd(:,2)&
1039 &   +deloc%rshift(3,s3)*rprimd(:,3)
1040    iprim=iprim+1
1041    call dang_d1(r1,r2,r3,bb)
1042    b_matrix(iprim,3*(i1-1)+1:3*i1) = b_matrix(iprim,3*(i1-1)+1:3*i1) + bb(:)
1043    call dang_d2(r1,r2,r3,bb)
1044    b_matrix(iprim,3*(i2-1)+1:3*i2) = b_matrix(iprim,3*(i2-1)+1:3*i2) + bb(:)
1045    call dang_d1(r3,r2,r1,bb)
1046    b_matrix(iprim,3*(i3-1)+1:3*i3) = b_matrix(iprim,3*(i3-1)+1:3*i3) + bb(:)
1047  end do
1048 
1049 !third: dihedral values
1050  do idihed=1,deloc%ndihed
1051    i1 = deloc%dihedrals(1,1,idihed)
1052    s1 = deloc%dihedrals(2,1,idihed)
1053    r1(:) = xcart(:,i1)+deloc%rshift(1,s1)*rprimd(:,1)&
1054 &   +deloc%rshift(2,s1)*rprimd(:,2)&
1055 &   +deloc%rshift(3,s1)*rprimd(:,3)
1056    i2 = deloc%dihedrals(1,2,idihed)
1057    s2 = deloc%dihedrals(2,2,idihed)
1058    r2(:) = xcart(:,i2)+deloc%rshift(1,s2)*rprimd(:,1)&
1059 &   +deloc%rshift(2,s2)*rprimd(:,2)&
1060 &   +deloc%rshift(3,s2)*rprimd(:,3)
1061    i3 = deloc%dihedrals(1,3,idihed)
1062    s3 = deloc%dihedrals(2,3,idihed)
1063    r3(:) = xcart(:,i3)+deloc%rshift(1,s3)*rprimd(:,1)&
1064 &   +deloc%rshift(2,s3)*rprimd(:,2)&
1065 &   +deloc%rshift(3,s3)*rprimd(:,3)
1066    i4 = deloc%dihedrals(1,4,idihed)
1067    s4 = deloc%dihedrals(2,4,idihed)
1068    r4(:) = xcart(:,i4)+deloc%rshift(1,s4)*rprimd(:,1)&
1069 &   +deloc%rshift(2,s4)*rprimd(:,2)&
1070 &   +deloc%rshift(3,s4)*rprimd(:,3)
1071 !  write(std_out,*) 'dihed ',idihed
1072 !  write(std_out,*) r1
1073 !  write(std_out,*) r2
1074 !  write(std_out,*) r3
1075 !  write(std_out,*) r4
1076 
1077    iprim=iprim+1
1078    call ddihedral_d1(r1,r2,r3,r4,bb)
1079    b_matrix(iprim,3*(i1-1)+1:3*i1) = b_matrix(iprim,3*(i1-1)+1:3*i1) + bb(:)
1080    call ddihedral_d2(r1,r2,r3,r4,bb)
1081    b_matrix(iprim,3*(i2-1)+1:3*i2) = b_matrix(iprim,3*(i2-1)+1:3*i2) + bb(:)
1082    call ddihedral_d2(r4,r3,r2,r1,bb)
1083    b_matrix(iprim,3*(i3-1)+1:3*i3) = b_matrix(iprim,3*(i3-1)+1:3*i3) + bb(:)
1084    call ddihedral_d1(r4,r3,r2,r1,bb)
1085    b_matrix(iprim,3*(i4-1)+1:3*i4) = b_matrix(iprim,3*(i4-1)+1:3*i4) + bb(:)
1086  end do
1087 
1088  do icart=1,deloc%ncart
1089    iprim=iprim+1
1090    b_matrix(iprim,3*(deloc%carts(2,icart)-1)+deloc%carts(1,icart)) = &
1091 &   b_matrix(iprim,3*(deloc%carts(2,icart)-1)+deloc%carts(1,icart)) + one
1092  end do
1093 
1094 !DEBUG
1095 ! write (200,*) 'calc_b_matrix : b_matrix = '
1096 ! do iprim=1,deloc%ninternal
1097 !   do i1=1, 3*natom
1098 !     write (200,'(E16.6,2x)',ADVANCE='NO') b_matrix(iprim,i1)
1099 !   end do
1100 !   write (200,*)
1101 ! end do
1102 !ENDDEBUG
1103 
1104 end subroutine calc_b_matrix

ABINIT/calc_btinv_matrix [ Functions ]

[ Top ] [ Functions ]

NAME

 calc_btinv_matrix

FUNCTION

INPUTS

OUTPUT

PARENTS

      xcart2deloc

NOTES

   bt_inv_matrix is inverse transpose of the delocalized
    coordinate B matrix. b_matrix is the primitive internal B matrix

CHILDREN

SOURCE

1671  subroutine calc_btinv_matrix(b_matrix,natom,ninternal,bt_inv_matrix,u_matrix)
1672 
1673 
1674 !This section has been created automatically by the script Abilint (TD).
1675 !Do not modify the following lines by hand.
1676 #undef ABI_FUNC
1677 #define ABI_FUNC 'calc_btinv_matrix'
1678 !End of the abilint section
1679 
1680  implicit none
1681 
1682 !Arguments ------------------------------------
1683  integer,intent(in) :: ninternal,natom
1684  real(dp),intent(in) :: b_matrix(ninternal,3*natom)
1685  real(dp),intent(out) :: bt_inv_matrix(3*(natom-1),3*natom)
1686  real(dp),intent(inout) :: u_matrix(ninternal,3*(natom-1))
1687 
1688 !Local variables ------------------------------------
1689 !scalars
1690  integer :: ii,info,lwork
1691 !arrays
1692  real(dp) :: f_eigs(3*natom),f_matrix(3*natom,3*natom)
1693  real(dp) :: s_matrix(3*natom,3*natom)
1694  real(dp) :: s_red(3*natom,3*(natom-1))
1695  real(dp) :: u_matrix_old(ninternal,3*(natom-1))
1696  real(dp),allocatable :: work(:)
1697 
1698 !******************************************************************
1699 
1700 !f matrix = B^{T} B
1701  call dgemm('T','N',3*natom,3*natom,ninternal,one,&
1702 & b_matrix,ninternal,b_matrix,ninternal,zero,f_matrix,3*natom)
1703 
1704  lwork = max(1,3*3*natom-1)
1705  ABI_ALLOCATE(work,(lwork))
1706  s_matrix(:,:) = f_matrix(:,:)
1707 
1708  call dsyev('V','L',3*natom,s_matrix,3*natom,f_eigs,work,lwork,info)
1709 
1710  ABI_DEALLOCATE(work)
1711 
1712  if (abs(f_eigs(1)) + abs(f_eigs(2)) + abs(f_eigs(3)) > tol10 ) then
1713    write(std_out,*) 'Error: 3 lowest eigenvalues are not zero'
1714    write(std_out,*) '  internal coordinates do NOT span the full degrees of freedom !'
1715    write(std_out,'(6E16.6)') f_eigs
1716    MSG_ERROR("Aborting now")
1717  end if
1718  if ( abs(f_eigs(4)) < tol10 ) then
1719    write(std_out,*) 'Error: fourth eigenvalue is zero'
1720    write(std_out,*) '  internal coordinates do NOT span the full degrees of freedom !'
1721    write(std_out,'(6E16.6)') f_eigs
1722    MSG_ERROR("Aborting now")
1723  end if
1724 
1725 !calculate U matrix from U = B * S_red * lambda^{-1/2}
1726  do ii=1,3*(natom-1)
1727    s_red(:,ii) = s_matrix(:,ii+3)/sqrt(f_eigs(ii+3))
1728  end do
1729 
1730  u_matrix_old(:,:) = u_matrix(:,:)
1731 
1732  call dgemm('N','N',ninternal,3*(natom-1),3*natom,one,&
1733 & b_matrix,ninternal,s_red,3*natom,zero,u_matrix,ninternal)
1734 
1735 
1736 !align eigenvectors, to preserve a form of continuity in convergences
1737 !!!! eigenvalues are no longer in increasing order!!! but only s_red is reordered
1738 !so that btinv is correct.
1739  call align_u_matrices(natom,ninternal,u_matrix,u_matrix_old,s_matrix,f_eigs)
1740 
1741 !calculate B_deloc^{-1} matrix for transformation of forces to deloc coord.
1742 !(B^{T}_deloc)^{-1} = (B_deloc B^{T}_deloc)^{-1} B_deloc = lambda^{-3/2} S^{T} F
1743 != ( S lambda^{3/2} )^{T} F
1744 
1745 !! DEFINITION
1746 !! real(dp),intent(out) :: bt_inv_matrix(3*(natom-1),3*natom)
1747 
1748 !even better: B_deloc^{-1} = lambda^{-1/2} S^{T}
1749  do ii=1,3*(natom-1)
1750 !  s_red(:,ii) = s_matrix(:,ii+3)*sqrt(f_eigs(ii+3))
1751    bt_inv_matrix(ii,:) = s_matrix(:,ii+3)/sqrt(f_eigs(ii+3))
1752  end do
1753 
1754 end subroutine calc_btinv_matrix

ABINIT/dang_d1 [ Functions ]

[ Top ] [ Functions ]

NAME

 dang_d1

FUNCTION

PARENTS

      calc_b_matrix

CHILDREN

      acrossb

SOURCE

1165 subroutine dang_d1(r1,r2,r3,bb)
1166 
1167 
1168 !This section has been created automatically by the script Abilint (TD).
1169 !Do not modify the following lines by hand.
1170 #undef ABI_FUNC
1171 #define ABI_FUNC 'dang_d1'
1172 !End of the abilint section
1173 
1174  implicit none
1175 
1176 !Arguments ------------------------------------
1177 !arrays
1178  real(dp),intent(in) :: r1(3),r2(3),r3(3)
1179  real(dp),intent(out) :: bb(3)
1180 
1181 !Local variables ------------------------------
1182 !scalars
1183  real(dp) :: cos_ang,n1,n1232,n2,tmp
1184 !arrays
1185  real(dp) :: cp1232(3),rpt(3),rpt12(3),rpt32(3)
1186 
1187 !************************************************************************
1188  n1=bond_length(r1,r2)
1189  n2=bond_length(r3,r2)
1190 
1191  rpt12(:) = r1(:)-r2(:)
1192  rpt32(:) = r3(:)-r2(:)
1193 
1194  cos_ang = (rpt12(1)*rpt32(1)+rpt12(2)*rpt32(2)+rpt12(3)*rpt32(3))/n1/n2
1195  if (cos_ang > one - epsilon(one)*two) then
1196    cos_ang = one
1197  else if(cos_ang < -one + epsilon(one)*two) then
1198    cos_ang = -one
1199  end if
1200 
1201  rpt(:) = rpt32(:)/n1/n2 - rpt12(:)*cos_ang/n1/n1
1202 
1203  tmp = sqrt(one-cos_ang**2)
1204  bb(:) = zero
1205  if (tmp > epsilon(one)) then
1206    bb(:) = rpt(:) * (-one)/tmp
1207  end if
1208 
1209 !TEST: version from MOLECULAR VIBRATIONS EB Wilson
1210  call acrossb(rpt12,rpt32,cp1232)
1211  n1232 = sqrt(cp1232(1)**2+cp1232(2)**2+cp1232(3)**2)
1212  rpt(:) = (cos_ang*rpt12(:)*n2/n1 - rpt32(:))/n1232
1213  if (abs(bb(1)-rpt(1))+abs(bb(2)-rpt(2))+abs(bb(3)-rpt(3)) > tol10) then
1214    write(std_out,*) 'Compare bb ang 1 : '
1215    write(std_out,*) bb(:), rpt(:), bb(:)-rpt(:)
1216  end if
1217  bb(:) = rpt(:)
1218 
1219 end subroutine dang_d1

ABINIT/dang_d2 [ Functions ]

[ Top ] [ Functions ]

NAME

 dang_d2

FUNCTION

PARENTS

      calc_b_matrix

CHILDREN

      acrossb

SOURCE

1238 subroutine dang_d2(r1,r2,r3,bb)
1239 
1240 
1241 !This section has been created automatically by the script Abilint (TD).
1242 !Do not modify the following lines by hand.
1243 #undef ABI_FUNC
1244 #define ABI_FUNC 'dang_d2'
1245 !End of the abilint section
1246 
1247  implicit none
1248 
1249 !Arguments ------------------------------------
1250 !arrays
1251  real(dp),intent(in) :: r1(3),r2(3),r3(3)
1252  real(dp),intent(out) :: bb(3)
1253 
1254 !Local variables ------------------------------
1255 !scalars
1256  real(dp) :: cos_ang,n1,n1232,n2,tmp
1257 !arrays
1258  real(dp) :: cp1232(3),rpt(3),rpt12(3),rpt32(3)
1259 
1260 !************************************************************************
1261  n1=bond_length(r1,r2)
1262  n2=bond_length(r3,r2)
1263 
1264  rpt12(:) = r1(:)-r2(:)
1265  rpt32(:) = r3(:)-r2(:)
1266 
1267  cos_ang = (rpt12(1)*rpt32(1)+rpt12(2)*rpt32(2)+rpt12(3)*rpt32(3))/n1/n2
1268  if (cos_ang > one - epsilon(one)*two) then
1269    cos_ang = one
1270  else if(cos_ang < -one + epsilon(one)*two) then
1271    cos_ang = -one
1272  end if
1273 
1274  rpt(:) = -rpt32(:)/n1/n2 - rpt12(:)/n1/n2 &
1275 & + rpt12(:)*cos_ang/n1/n1 + rpt32(:)*cos_ang/n2/n2
1276 
1277  tmp = sqrt(one-cos_ang**2)
1278  bb(:) = zero
1279  if (tmp > tol12) then
1280    bb(:) = rpt(:) * (-one)/tmp
1281  end if
1282 
1283 !TEST: version from MOLECULAR VIBRATIONS EB Wilson
1284  call acrossb(rpt12,rpt32,cp1232)
1285  n1232 = sqrt(cp1232(1)**2+cp1232(2)**2+cp1232(3)**2)
1286  rpt(:) = ((n1-n2*cos_ang)*rpt12(:)/n1 + (n2-n1*cos_ang)*rpt32(:)/n2) / n1232
1287  if (abs(bb(1)-rpt(1))+abs(bb(2)-rpt(2))+abs(bb(3)-rpt(3))  > tol10) then
1288    write(std_out,*) 'Compare bb ang 2 : '
1289    write(std_out,*) bb(:), rpt(:), bb(:)-rpt(:)
1290  end if
1291  bb(:) = rpt(:)
1292 
1293 end subroutine dang_d2

ABINIT/dbond_length_d1 [ Functions ]

[ Top ] [ Functions ]

NAME

 dbond_length_d1

FUNCTION

PARENTS

      calc_b_matrix

CHILDREN

      acrossb

SOURCE

1122 subroutine dbond_length_d1(r1,r2,bb)
1123 
1124 
1125 !This section has been created automatically by the script Abilint (TD).
1126 !Do not modify the following lines by hand.
1127 #undef ABI_FUNC
1128 #define ABI_FUNC 'dbond_length_d1'
1129 !End of the abilint section
1130 
1131  implicit none
1132 
1133 !Arguments ------------------------------------
1134 !arrays
1135  real(dp),intent(in) :: r1(3),r2(3)
1136  real(dp),intent(out) :: bb(3)
1137 
1138 !Local variables ------------------------------
1139 !arrays
1140  real(dp) :: rpt(3)
1141 
1142 !************************************************************************
1143  rpt(:) = r1(:)-r2(:)
1144  bb(:) = rpt(:)/bond_length(r1,r2)
1145 
1146 end subroutine dbond_length_d1

ABINIT/ddihedral_d1 [ Functions ]

[ Top ] [ Functions ]

NAME

 ddihedral_d1

FUNCTION

PARENTS

      calc_b_matrix

CHILDREN

      acrossb

SOURCE

1311 subroutine ddihedral_d1(r1,r2,r3,r4,bb)
1312 
1313 
1314 !This section has been created automatically by the script Abilint (TD).
1315 !Do not modify the following lines by hand.
1316 #undef ABI_FUNC
1317 #define ABI_FUNC 'ddihedral_d1'
1318 !End of the abilint section
1319 
1320  implicit none
1321 
1322 !Arguments ------------------------------------
1323 !arrays
1324  real(dp),intent(in) :: r1(3),r2(3),r3(3),r4(3)
1325  real(dp),intent(out) :: bb(3)
1326 
1327 !Local variables ------------------------------------
1328 !scalars
1329  real(dp) :: cos_dihedral,dih_sign,n1,n2,n23,sin_dihedral,tmp
1330 !arrays
1331  real(dp) :: cp1232(3),cp32_1232(3),cp32_3432(3),cp3432(3),cpcp(3),rpt(3)
1332  real(dp) :: rpt12(3),rpt32(3),rpt34(3)
1333 
1334 !******************************************************************
1335  rpt12(:) = r1(:)-r2(:)
1336  rpt32(:) = r3(:)-r2(:)
1337  rpt34(:) = r3(:)-r4(:)
1338 
1339  call acrossb(rpt12,rpt32,cp1232)
1340  call acrossb(rpt34,rpt32,cp3432)
1341 
1342 !DEBUG
1343 !write(std_out,*) ' cos_dihedral : cp1232 = ', cp1232
1344 !write(std_out,*) ' cos_dihedral : cp3432 = ', cp3432
1345 !ENDDEBUG
1346 
1347  n1 = sqrt(cp1232(1)**2+cp1232(2)**2+cp1232(3)**2)
1348  n2 = sqrt(cp3432(1)**2+cp3432(2)**2+cp3432(3)**2)
1349 
1350  cos_dihedral = (cp1232(1)*cp3432(1)+cp1232(2)*cp3432(2)+cp1232(3)*cp3432(3))/n1/n2
1351  if (cos_dihedral > one - epsilon(one)*two) then
1352    cos_dihedral = one
1353  else if(cos_dihedral < -one + epsilon(one)*two) then
1354    cos_dihedral = -one
1355  end if
1356 !we use complementary of standard angle, so
1357 !cos_dihedral = -cos_dihedral
1358 
1359  call acrossb(cp1232,cp3432,cpcp)
1360  cpcp(:) = cpcp(:)/n1/n2
1361 !we use complementary of standard angle, but sin is invariant
1362  sin_dihedral = -(cpcp(1)*rpt32(1)+cpcp(2)*rpt32(2)+cpcp(3)*rpt32(3))&
1363 & /sqrt(rpt32(1)**2+rpt32(2)**2+rpt32(3)**2)
1364  dih_sign = one
1365  if (sin_dihedral < -epsilon(one)) then
1366    dih_sign = -one
1367  end if
1368 
1369 !DEBUG
1370 !write(std_out,'(a,3E16.6)') 'ddihedral_d1 : cos abs(sin) dih_sign= ',&
1371 !&    cos_dihedral,sin_dihedral,dih_sign
1372 !ENDDEBUG
1373 
1374 !ddihedral_d1 = dih_sign* acos(cos_dihedral)
1375  call acrossb(rpt32,cp1232,cp32_1232)
1376  call acrossb(rpt32,cp3432,cp32_3432)
1377 
1378  rpt(:) = cp32_3432(:)/n1/n2 - cp32_1232(:)/n1/n1 * cos_dihedral
1379  bb(:) = zero
1380 
1381 !DEBUG
1382 !write(std_out,*) 'ddihedral_d1 cp1232 cp3432 = ',cp1232,cp3432,rpt32
1383 !write(std_out,*) 'ddihedral_d1 cp32_1232 cp32_3432 = ',cp32_1232,cp32_3432,cos_dihedral,n1,n2
1384 !write(std_out,*) 'ddihedral_d1 rpt = ',rpt
1385 !ENDDEBUG
1386 
1387  tmp = sqrt(one-cos_dihedral**2)
1388  if (tmp > tol12) then
1389 !  we use complementary of standard angle, so cosine in acos has - sign,
1390 !  and it appears for the derivative
1391    bb(:) = -dih_sign * rpt(:) * (-one) / tmp
1392  else
1393    bb(:) = dih_sign * cp32_3432(:) / n1 / n2 / &
1394 &   sqrt(cp32_3432(1)**2+cp32_3432(2)**2+cp32_3432(3)**2)
1395  end if
1396 
1397 !TEST: version from MOLECULAR VIBRATIONS EB Wilson
1398 
1399  n23 = sqrt(rpt32(1)*rpt32(1)+rpt32(2)*rpt32(2)+rpt32(3)*rpt32(3))
1400  rpt(:) = cp1232(:)*n23/n1/n1
1401 !if (abs(bb(1)-rpt(1))+abs(bb(2)-rpt(2))+abs(bb(3)-rpt(3))  > tol10) then
1402 !write(std_out,*) 'Compare bb1 : '
1403 !write(std_out,*) bb(:), rpt(:), bb(:)-rpt(:)
1404 !end if
1405  bb(:) = rpt(:)
1406 
1407 end subroutine ddihedral_d1

ABINIT/ddihedral_d2 [ Functions ]

[ Top ] [ Functions ]

NAME

 ddihedral_d2

FUNCTION

PARENTS

      calc_b_matrix

CHILDREN

      acrossb

SOURCE

1425 subroutine ddihedral_d2(r1,r2,r3,r4,bb)
1426 
1427 
1428 !This section has been created automatically by the script Abilint (TD).
1429 !Do not modify the following lines by hand.
1430 #undef ABI_FUNC
1431 #define ABI_FUNC 'ddihedral_d2'
1432 !End of the abilint section
1433 
1434  implicit none
1435 
1436 !Arguments ------------------------------------
1437 !arrays
1438  real(dp),intent(in) :: r1(3),r2(3),r3(3),r4(3)
1439  real(dp),intent(out) :: bb(3)
1440 
1441 !Local variables
1442 !scalars
1443  real(dp) :: cos_dihedral,dih_sign,n1,n2,n23,sin_dihedral,sp1232,sp3432,tmp
1444 !arrays
1445  real(dp) :: cp1232(3),cp1232_12(3),cp1232_34(3),cp32_1232(3),cp32_3432(3)
1446  real(dp) :: cp3432(3),cp3432_12(3),cp3432_34(3),cpcp(3),rpt(3),rpt12(3)
1447  real(dp) :: rpt32(3),rpt34(3)
1448 
1449 ! *************************************************************************
1450  rpt12(:) = r1(:)-r2(:)
1451  rpt32(:) = r3(:)-r2(:)
1452  rpt34(:) = r3(:)-r4(:)
1453 
1454  call acrossb(rpt12,rpt32,cp1232)
1455  call acrossb(rpt34,rpt32,cp3432)
1456 
1457 !DEBUG
1458 !write(std_out,*) ' cos_dihedral : cp1232 = ', cp1232
1459 !write(std_out,*) ' cos_dihedral : cp3432 = ', cp3432
1460 !ENDDEBUG
1461 
1462  n1 = sqrt(cp1232(1)**2+cp1232(2)**2+cp1232(3)**2)
1463  n2 = sqrt(cp3432(1)**2+cp3432(2)**2+cp3432(3)**2)
1464 
1465  cos_dihedral = (cp1232(1)*cp3432(1)+cp1232(2)*cp3432(2)+cp1232(3)*cp3432(3))/n1/n2
1466  if (cos_dihedral > one - epsilon(one)*two) then
1467    cos_dihedral = one
1468  else if(cos_dihedral < -one + epsilon(one)*two) then
1469    cos_dihedral = -one
1470  end if
1471 !we use complementary of standard angle, so
1472 !cos_dihedral = -cos_dihedral
1473 
1474  call acrossb(cp1232,cp3432,cpcp)
1475  cpcp(:) = cpcp(:)/n1/n2
1476 !we use complementary of standard angle, but sin is invariant
1477  sin_dihedral = -(cpcp(1)*rpt32(1)+cpcp(2)*rpt32(2)+cpcp(3)*rpt32(3))&
1478 & /sqrt(rpt32(1)**2+rpt32(2)**2+rpt32(3)**2)
1479  dih_sign = one
1480  if (sin_dihedral <  -tol12) then
1481    dih_sign = -one
1482  end if
1483 
1484 !DEBUG
1485 !write(std_out,'(a,3E16.6)') 'ddihedral_d2 : cos abs(sin) dih_sign= ',&
1486 !&    cos_dihedral,sin_dihedral,dih_sign
1487 !ENDDEBUG
1488 
1489 !ddihedral_d2 = dih_sign* acos(cos_dihedral)
1490  call acrossb(rpt32,cp3432,cp32_3432)
1491  call acrossb(cp3432,rpt12,cp3432_12)
1492  call acrossb(cp1232,rpt34,cp1232_34)
1493 
1494  call acrossb(rpt32,cp1232,cp32_1232)
1495  call acrossb(cp1232,rpt12,cp1232_12)
1496  call acrossb(cp3432,rpt34,cp3432_34)
1497 
1498  rpt(:) = -(cp32_3432(:) + cp3432_12(:) + cp1232_34(:))/n1/n2 &
1499 & +cos_dihedral*(cp32_1232(:)/n1/n1 + cp1232_12(:)/n1/n1 + cp3432_34(:)/n2/n2)
1500  bb(:) = zero
1501  tmp = sqrt(one-cos_dihedral**2)
1502  if (tmp > tol12) then
1503 !  we use complementary of standard angle, so cosine in acos has - sign,
1504 !  and it appears for derivative
1505    bb(:) = -dih_sign * rpt(:) * (-one) / tmp
1506  else
1507    bb(:) = dih_sign * cos_dihedral * &
1508 &   ( cp32_1232(:)/n1/n1/sqrt(cp32_1232(1)**2+cp32_1232(2)**2+cp32_1232(3)**2) &
1509 &   +cp1232_12(:)/n1/n1/sqrt(cp1232_12(1)**2+cp1232_12(2)**2+cp1232_12(3)**2) &
1510 &   +cp3432_34(:)/n2/n2/sqrt(cp3432_34(1)**2+cp3432_34(2)**2+cp3432_34(3)**2) )
1511  end if
1512 
1513 !TEST: version from MOLECULAR VIBRATIONS EB Wilson p. 61
1514  n23 = sqrt(rpt32(1)*rpt32(1)+rpt32(2)*rpt32(2)+rpt32(3)*rpt32(3))
1515  sp1232 = rpt12(1)*rpt32(1)+rpt12(2)*rpt32(2)+rpt12(3)*rpt32(3)
1516  sp3432 = rpt34(1)*rpt32(1)+rpt34(2)*rpt32(2)+rpt34(3)*rpt32(3)
1517 
1518  rpt(:) = -cp1232(:)*(n23-sp1232/n23)/n1/n1 - cp3432(:)*sp3432/n23/n2/n2
1519 !if (abs(bb(1)-rpt(1))+abs(bb(2)-rpt(2))+abs(bb(3)-rpt(3))  > tol10) then
1520 !write(std_out,*) 'Compare bb2 : '
1521 !write(std_out,*) bb(:), rpt(:), bb(:)-rpt(:)
1522 !write(std_out,*) -cp1232(:)*(n23-sp1232/n23)/n1/n1, -cp3432(:)*sp3432/n23/n2/n2
1523 !end if
1524  bb(:) = rpt(:)
1525 
1526 end subroutine ddihedral_d2

ABINIT/deloc2xcart [ Functions ]

[ Top ] [ 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

PARENTS

      pred_delocint

CHILDREN

      dgemv,wrtout,xcart2deloc

SOURCE

705 subroutine deloc2xcart(deloc,natom,rprimd,xcart,deloc_int,btinv,u_matrix)
706 
707 
708 !This section has been created automatically by the script Abilint (TD).
709 !Do not modify the following lines by hand.
710 #undef ABI_FUNC
711 #define ABI_FUNC 'deloc2xcart'
712 !End of the abilint section
713 
714  implicit none
715 
716 !Arguments ------------------------------------
717 !scalars
718  integer,intent(in) :: natom
719  type(delocint),intent(in) :: deloc
720 !arrays
721  real(dp),intent(in) :: deloc_int(3*(natom-1)),rprimd(3,3)
722  real(dp),intent(inout) :: u_matrix(deloc%ninternal,3*(natom-1))
723  real(dp),intent(inout) :: xcart(3,natom)
724  real(dp),intent(out) :: btinv(3*(natom-1),3*natom)
725 
726 !Local variables-------------------------------
727 !scalars
728  integer :: iiter,iprim,niter
729  integer :: ii
730  real(dp) :: minmix, maxmix
731  real(dp) :: mix,tot_diff, toldeloc
732  real(dp) :: lntoldeloc
733  logical  :: DEBUG=.FALSE.
734 !arrays
735  real(dp) :: btinv_tmp(3*(natom-1),3*natom)
736  real(dp) :: cgrad(3*natom),cgrad_old(3*natom)
737  real(dp) :: deloc_int_now(3*(natom-1)),prim_int(deloc%ninternal)
738  real(dp) :: tmpxcart(3*natom)
739  real(dp) :: xdeloc_diff(3*(natom-1))
740 
741  character(len=500) :: message
742 
743 ! ******************************************************************
744 
745  if (DEBUG) then
746    write(ab_out,*) 'ENTERING DELOC2XCART'
747 
748    write (message,*) 'BONDS=',deloc%nbond
749    call wrtout(ab_out,message,'COLL')
750    do ii = 1, deloc%nbond
751      write (message,*) ii, deloc%bonds(:,:,ii)
752      call wrtout(ab_out,message,'COLL')
753    end do
754 
755    write (message,*) 'ANGS=',deloc%nang
756    call wrtout(ab_out,message,'COLL')
757    do ii = 1, deloc%nang
758      write (message,*) ii, deloc%angs(:,:,ii)
759      call wrtout(ab_out,message,'COLL')
760    end do
761 
762    write (message,*) 'DIHEDRALS=',deloc%ndihed
763    call wrtout(ab_out,message,'COLL')
764    do ii = 1, deloc%ndihed
765      write (message,*) ii, deloc%dihedrals(:,:,ii)
766      call wrtout(ab_out,message,'COLL')
767    end do
768 
769    write (message,*) 'CARTS=',deloc%ncart
770    call wrtout(ab_out,message,'COLL')
771    do ii = 1, deloc%ncart
772      write (message,*) ii, deloc%carts(:,ii)
773      call wrtout(ab_out,message,'COLL')
774    end do
775 
776    write (ab_out,*) 'xcart (input)'
777    do ii=1,natom
778      write (ab_out,*) xcart(:,ii)
779    end do
780 
781  end if
782 
783  niter = 200
784  tmpxcart = reshape(xcart,(/3*natom/))
785 
786  cgrad_old(:) = zero
787  cgrad(:) = zero
788  maxmix = 0.9_dp
789  minmix = 0.2_dp
790  toldeloc = tol10
791  lntoldeloc = log(toldeloc)
792 
793  do iiter=1,niter
794    if (iiter==1) then
795      mix= minmix
796    else
797      mix = minmix + (maxmix-minmix)*(log(tot_diff)-lntoldeloc) / lntoldeloc
798    end if
799    if (mix < minmix) mix = minmix
800    if (mix > maxmix) mix = maxmix
801 
802    tmpxcart(:) = tmpxcart(:) + mix*cgrad(:)
803    xcart = reshape(tmpxcart,(/3,natom/))
804    call xcart2deloc(deloc,natom,rprimd,xcart,&
805 &   btinv_tmp,u_matrix,deloc_int_now,prim_int)
806 !  update the BT^{-1} matrix?
807    btinv(:,:) = btinv_tmp(:,:)
808 
809    xdeloc_diff(:) = deloc_int(:) - deloc_int_now(:)
810 
811    tot_diff = sum(abs(xdeloc_diff))
812    if (tot_diff < toldeloc) exit
813 
814    cgrad_old(:) = cgrad(:)
815 
816 !  gradient vector = btinv^{T} * xdeloc_diff
817    call dgemv('T',3*(natom-1),3*natom,one,&
818 &   btinv,3*(natom-1),xdeloc_diff,1,zero,cgrad,1)
819  end do
820 !end iiter do
821 
822  call xcart2deloc(deloc,natom,rprimd,xcart,&
823 & btinv,u_matrix,deloc_int_now,prim_int)
824  write (message,'(3a)') 'delocalized internals, after convergence of xcart = ', ch10
825  call wrtout(std_out,message,'COLL')
826  do ii = 1, 3*(natom-1)
827    write (message,'(I6,E20.10,2x)') ii, deloc_int_now(ii)
828    call wrtout(std_out,message,'COLL')
829  end do
830 
831  xdeloc_diff(:) = deloc_int(:) - deloc_int_now(:)
832 
833  write (message,'(a)') 'Primitive internal coordinate values:'
834  call wrtout(std_out,message,'COLL')
835  do iprim = 1, deloc%nbond
836    write (message,'(i6,E20.10)') iprim, prim_int(iprim)
837    call wrtout(std_out,message,'COLL')
838  end do
839  do iprim = deloc%nbond+1, deloc%nbond+deloc%nang+deloc%ndihed
840    write (message,'(i6,2E20.10)') iprim, prim_int(iprim), prim_int(iprim)/pi*180.0_dp
841    call wrtout(std_out,message,'COLL')
842  end do
843  do iprim = deloc%nbond+deloc%nang+deloc%ndihed+1, deloc%ninternal
844    write (message,'(i6,E20.10)') iprim, prim_int(iprim)
845    call wrtout(std_out,message,'COLL')
846  end do
847 
848  if (iiter == niter+1) then
849    write (message,'(a,i6,a,E20.10)') 'deloc2xcart : Error, xcart not converged in ', niter, 'iterations ', tot_diff
850    MSG_ERROR(message)
851  end if
852 
853  if(DEBUG)then
854    write (ab_out,*) 'xcart (output)'
855    do ii=1,natom
856      write (ab_out,*) xcart(:,ii)
857    end do
858    write(ab_out,*) 'EXITING DELOC2XCART'
859  end if
860 
861 end subroutine deloc2xcart

ABINIT/fred2fdeloc [ Functions ]

[ Top ] [ Functions ]

NAME

 fred2fdeloc

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_force(3*(natom-1))=delocalized forces from reduced coordinate ones
 fred(3,natom)=delocalized forces in reduced coordinates

PARENTS

      pred_delocint,xfh_recover_deloc

CHILDREN

      dgemm,dgemv,wrtout

SOURCE

888 subroutine fred2fdeloc(btinv,deloc_force,fred,natom,gprimd)
889 
890 
891 !This section has been created automatically by the script Abilint (TD).
892 !Do not modify the following lines by hand.
893 #undef ABI_FUNC
894 #define ABI_FUNC 'fred2fdeloc'
895 !End of the abilint section
896 
897  implicit none
898 
899 !Arguments ------------------------------------
900 !scalars
901  integer, intent(in) :: natom
902 !arrays
903  real(dp),intent(in) :: btinv(3*(natom-1),3*natom),gprimd(3,3),fred(3,natom)
904  real(dp),intent(out) :: deloc_force(3*(natom-1))
905 
906 !Local variables-------------------------------
907  integer :: ii
908 !arrays
909  real(dp) :: fcart(3,natom)
910  character(len=500) :: message
911 
912 ! ******************************************************************
913 
914 !make cartesian forces
915 
916  call dgemm('N','N',3,natom,3,one,&
917 & gprimd,3,fred,3,zero,fcart,3)
918 
919 !turn cartesian to delocalized forces
920  call dgemv('N',3*(natom-1),3*natom,one,&
921 & btinv,3*(natom-1),fcart,1,zero,deloc_force,1)
922 
923  write (message,'(a)') 'fred2fdeloc : deloc_force = '
924  call wrtout(std_out,message,'COLL')
925 
926  do ii = 1, 3*(natom-1)
927    write (message,'(I6,E16.6)') ii, deloc_force(ii)
928    call wrtout(std_out,message,'COLL')
929  end do
930 
931 end subroutine fred2fdeloc

ABINIT/m_pred_delocint [ Modules ]

[ Top ] [ Modules ]

NAME

 m_pred_delocint

FUNCTION

COPYRIGHT

  Copyright (C) 1998-2018 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 .

PARENTS

CHILDREN

SOURCE

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 #include "abi_common.h"
25 
26 module m_pred_delocint
27 
28  use defs_basis
29  use m_errors
30  use m_abicore
31  use m_abimover
32  use m_abihist
33  use m_xfpack
34  use m_linalg_interfaces
35 
36  use m_geometry,   only : fcart2fred, xcart2xred, xred2xcart, metric, acrossb
37  use m_bfgs,       only : hessinit, hessupdt, brdene
38  use m_results_gs, only : results_gs_type
39 
40  implicit none
41 
42  private

ABINIT/pred_delocint [ Functions ]

[ Top ] [ 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

PARENTS

      mover

CHILDREN

      brdene,calc_prim_int,deloc2xcart,fcart2fred,fred2fdeloc,hessupdt
      hist2var,make_prim_internals,metric,var2hist,wrtout,xcart2deloc
      xcart2xred,xfh_recover_deloc,xfpack_f2vout,xfpack_vin2x,xfpack_x2vin
      xred2xcart

SOURCE

103 subroutine pred_delocint(ab_mover,ab_xfh,deloc,forstr,hist,ionmov,itime,zDEBUG,iexit)
104 
105 
106 !This section has been created automatically by the script Abilint (TD).
107 !Do not modify the following lines by hand.
108 #undef ABI_FUNC
109 #define ABI_FUNC 'pred_delocint'
110 !End of the abilint section
111 
112  implicit none
113 
114 !Arguments ------------------------------------
115 !scalars
116  type(abimover),intent(in)       :: ab_mover
117  type(ab_xfh_type),intent(inout)    :: ab_xfh
118  type(abihist),intent(inout) :: hist
119  type(abiforstr),intent(in) :: forstr
120  type(delocint),intent(inout) :: deloc
121  integer,intent(in) :: itime
122  integer,intent(in) :: ionmov
123  integer,intent(in) :: iexit
124  logical,intent(in) :: zDEBUG
125 
126 !Local variables-------------------------------
127 !scalars
128  integer  :: ndim,cycl_main
129  integer  :: ihist_prev,ii,jj,kk
130  real(dp),save :: ucvol0
131  real(dp) :: ucvol
132  real(dp) :: etotal,etotal_prev
133  logical  :: DEBUG=.TRUE.
134  integer,save :: icenter,irshift ! DELOCINT indexes
135  integer,save :: nshell,ndeloc ! DELOCINT number of
136  character(len=500) :: message
137 
138 !arrays
139  real(dp),allocatable,save :: hessin(:,:),vin(:),vin_prev(:)
140  real(dp),allocatable,save :: vout(:),vout_prev(:)
141  real(dp),save :: acell0(3) ! Initial acell
142  real(dp),save :: rprimd0(3,3) ! Initial rprimd
143  real(dp),allocatable :: prim_int(:)
144  real(dp),allocatable,save :: u_matrix(:,:) ! DELOCINT this may need to be added to type inside ab_mover
145  real(dp) :: acell(3)
146  real(dp) :: rprimd(3,3)
147  real(dp) :: gprimd(3,3)
148  real(dp) :: gmet(3,3)
149  real(dp) :: rmet(3,3)
150  real(dp) :: residual(3,ab_mover%natom)
151 !real(dp) :: residual_corrected(3,ab_mover%natom)
152  real(dp) :: xred(3,ab_mover%natom),xcart(3,ab_mover%natom)
153  real(dp) :: strten(6)
154  real(dp) :: deloc_force(3*(ab_mover%natom-1))
155  real(dp) :: deloc_int(3*(ab_mover%natom-1))
156  real(dp) :: bt_inv_matrix(3*(ab_mover%natom-1),3*ab_mover%natom)
157 
158 !***************************************************************************
159 !Beginning of executable session
160 !***************************************************************************
161 
162  if(iexit/=0)then
163    if (allocated(vin))        then
164      ABI_DEALLOCATE(vin)
165    end if
166    if (allocated(vout))       then
167      ABI_DEALLOCATE(vout)
168    end if
169    if (allocated(vin_prev))   then
170      ABI_DEALLOCATE(vin_prev)
171    end if
172    if (allocated(vout_prev))  then
173      ABI_DEALLOCATE(vout_prev)
174    end if
175    if (allocated(hessin))     then
176      ABI_DEALLOCATE(hessin)
177    end if
178    if (allocated(u_matrix))     then
179      ABI_DEALLOCATE(u_matrix)
180    end if
181    return
182  end if
183 
184 !write(std_out,*) 'delocint 01'
185 !##########################################################
186 !### 01. Debugging and Verbose
187 
188  if(DEBUG)then
189    write(std_out,'(a,3a,38a,39a)') ch10,('-',kk=1,3),&
190 &   'Debugging and Verbose for pred_deloint',('-',kk=1,39)
191    write(std_out,*) 'ionmov: ',ionmov
192    write(std_out,*) 'itime:  ',itime
193  end if
194 
195 !write(std_out,*) 'delocint 02'
196 !##########################################################
197 !### 02. Compute the dimension of vectors (ndim)
198 
199 !With internal we have 1 coordinate less
200  ndeloc = 3*(ab_mover%natom-1)
201  ndim=ndeloc
202  deloc_int(:)=zero
203  deloc_force(:)=zero
204  if(ab_mover%optcell==1 .or.&
205 & ab_mover%optcell==4 .or.&
206 & ab_mover%optcell==5 .or.&
207 & ab_mover%optcell==6) ndim=ndim+1
208  if(ab_mover%optcell==2 .or.&
209 & ab_mover%optcell==3) ndim=ndim+6
210  if(ab_mover%optcell==7 .or.&
211 & ab_mover%optcell==8 .or.&
212 & ab_mover%optcell==9) ndim=ndim+3
213 
214  if(DEBUG) write(std_out,*) 'Dimension of vin, vout and hessian (ndim): ',ndim
215 
216 !write(std_out,*) 'delocint 03'
217 !##########################################################
218 !### 03. Allocate the vectors vin, vout and hessian matrix
219 
220 !Notice thqt vin, vout, etc could be allocated
221 !From a previous dataset with a different ndim
222  if(itime==1)then
223    if (allocated(vin))        then
224      ABI_DEALLOCATE(vin)
225    end if
226    if (allocated(vout))       then
227      ABI_DEALLOCATE(vout)
228    end if
229    if (allocated(vin_prev))   then
230      ABI_DEALLOCATE(vin_prev)
231    end if
232    if (allocated(vout_prev))  then
233      ABI_DEALLOCATE(vout_prev)
234    end if
235    if (allocated(hessin))     then
236      ABI_DEALLOCATE(hessin)
237    end if
238    ABI_ALLOCATE(vin,(ndim))
239    ABI_ALLOCATE(vout,(ndim))
240    ABI_ALLOCATE(vin_prev,(ndim))
241    ABI_ALLOCATE(vout_prev,(ndim))
242    ABI_ALLOCATE(hessin,(ndim,ndim))
243 
244  end if
245 
246 
247 !write(std_out,*) 'delocint 04'
248 !##########################################################
249 !### 04. Obtain the present values from the history
250 
251  call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
252  call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
253 
254  strten(:)=hist%strten(:,hist%ihist)
255  etotal   =hist%etot(hist%ihist)
256 
257 !Fill the residual with forces (No preconditioning)
258 !Or the preconditioned forces
259  if (ab_mover%goprecon==0)then
260    call fcart2fred(hist%fcart(:,:,hist%ihist),residual,rprimd,ab_mover%natom)
261  else
262    residual(:,:)= forstr%fred(:,:)
263  end if
264 
265  if(zDEBUG)then
266    write (std_out,*) 'residual:'
267    do kk=1,ab_mover%natom
268      write (std_out,*) residual(:,kk)
269    end do
270    write (std_out,*) 'strten:'
271    write (std_out,*) strten(1:3),ch10,strten(4:6)
272    write (std_out,*) 'etotal:'
273    write (std_out,*) etotal
274  end if
275 
276  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
277 
278 !Save initial values
279  if (itime==1)then
280    acell0(:)=acell(:)
281    rprimd0(:,:)=rprimd(:,:)
282    ucvol0=ucvol
283  end if
284 
285 !DEBUG (UCVOL)
286  if(DEBUG)then
287    write(std_out,*) 'Volume of cell (ucvol):',ucvol
288  end if
289 
290 !Get rid of mean force on whole unit cell, but only if no
291 !generalized constraints are in effect
292 !  residual_corrected(:,:)=residual(:,:)
293 !  if(ab_mover%nconeq==0)then
294 !    do ii=1,3
295 !      if (ii/=3.or.ab_mover%jellslab==0) then
296 !        favg=sum(residual_corrected(ii,:))/dble(ab_mover%natom)
297 !        residual_corrected(ii,:)=residual_corrected(ii,:)-favg
298 !      end if
299 !    end do
300 !  end if
301 
302 !write(std_out,*) 'delocint 05'
303 !##########################################################
304 !### 05. Compute internals for first time
305 
306  if (itime==1)then
307    call make_prim_internals(deloc,ab_mover%natom,&
308 &   ab_mover%ntypat,rprimd,ab_mover%typat,xcart,ab_mover%znucl)
309 
310    ABI_ALLOCATE(prim_int,(deloc%ninternal))
311 
312    if(DEBUG)then
313      write (message,'(a,i6)') 'Number of primitive internal coordinates (ninternal): ',deloc%ninternal
314      call wrtout(std_out,  message,'COLL')
315    end if
316 
317    if (allocated(u_matrix))  then
318      ABI_DEALLOCATE(u_matrix)
319    end if
320    ABI_ALLOCATE(u_matrix,(deloc%ninternal,ndeloc))
321 
322    call calc_prim_int(deloc,ab_mover%natom,rprimd,xcart,prim_int)
323 
324    if(DEBUG)then
325      write (message,'(a)') 'Primitive internal coordinate values:'
326      call wrtout(std_out,  message,'COLL')
327      write (message,'(a)') ' Bonds:'
328      call wrtout(std_out,  message,'COLL')
329      do ii = 1, deloc%nbond
330        write (message,'(i6,E20.10)') ii, prim_int(ii)
331        call wrtout(std_out,  message,'COLL')
332      end do
333 
334      write (message,'(a)') ' Angles:'
335      call wrtout(std_out,  message,'COLL')
336      do ii = deloc%nbond+1, deloc%nbond+deloc%nang
337        write (message,'(i6,2(E20.10,2x))') ii, prim_int(ii), prim_int(ii)/pi*180.0_dp
338        call wrtout(std_out,  message,'COLL')
339      end do
340 
341      write (message,'(a)') ' Dihedrals:'
342      call wrtout(std_out,  message,'COLL')
343      do ii = deloc%nbond+deloc%nang+1, deloc%nbond+deloc%nang+deloc%ndihed
344        write (message,'(i6,2(E20.10,2x))') ii, prim_int(ii), prim_int(ii)/pi*180.0_dp
345        call wrtout(std_out,  message,'COLL')
346      end do
347 
348      write (message,'(a)') ' Cartesian auxiliary coordinates for constraints:'
349      call wrtout(std_out,  message,'COLL')
350      do ii = deloc%nbond+deloc%nang+deloc%ndihed+1, deloc%ninternal
351        write (message,'(i6,E20.10)') ii, prim_int(ii)
352        call wrtout(std_out,  message,'COLL')
353      end do
354    end if
355 
356    ABI_DEALLOCATE(prim_int)
357 
358 !  equal weight on all internal coordinates as a starting point.
359    u_matrix(:,:) = one / dble (ndeloc)
360 
361 !  Zero the arrays before first use
362    deloc_force(:) = zero
363 
364  end if
365 
366  ABI_ALLOCATE(prim_int,(deloc%ninternal))
367 
368 !write(std_out,*) 'delocint 06'
369 !##########################################################
370 !### 06. Compute delocalized coordinates and forces
371 
372 !xcart ---> deloc_int
373 
374 !Convert positions to delocalized coordinates for next step
375  call xcart2deloc(deloc,ab_mover%natom,rprimd,xcart,&
376 & bt_inv_matrix,u_matrix,deloc_int,prim_int)
377 
378 !fred ---> deloc_force
379 
380 !Convert forces to delocalized coordinates for next step
381  call fred2fdeloc(bt_inv_matrix,deloc_force,residual,ab_mover%natom,gprimd)
382 
383 !write(std_out,*) 'delocint 07'
384 !##########################################################
385 !### 07. Fill the vectors vin and vout
386 
387 !DEBUG deloc_int and deloc_force before pack
388  if(DEBUG)then
389    write (std_out,*) 'Delocalized internals and forces (ndeloc):',ndeloc
390    write(std_out,*) 'deloc_int'
391    do ii=1,ndeloc,3
392      if (ii+2<=ndeloc)then
393        write(std_out,*) ii,deloc_int(ii:ii+2)
394      else
395        write(std_out,*) ii,deloc_int(ii:ndeloc)
396      end if
397    end do
398    write(std_out,*) 'deloc_force'
399    do ii=1,ndeloc,3
400      if (ii+2<=ndeloc)then
401        write(std_out,*) ii,deloc_force(ii:ii+2)
402      else
403        write(std_out,*) ii,deloc_force(ii:ndeloc)
404      end if
405    end do
406  end if
407 
408 !DELOCINT
409 !Instead of fred_corrected we use deloc_force
410 !Instead of xred e use deloc_int
411 !
412 !Initialize input vectors : first vin, then vout
413 !The values of vin from the previous iteration
414 !should be the same
415  call xfpack_x2vin(acell, acell0, ab_mover%natom-1, ndim,&
416 & ab_mover%nsym, ab_mover%optcell, rprimd, rprimd0,&
417 & ab_mover%symrel, ucvol, ucvol0, vin, deloc_int)
418 !end if
419 
420  call xfpack_f2vout(deloc_force, ab_mover%natom-1, ndim,&
421 & ab_mover%optcell, ab_mover%strtarget, strten, ucvol, &
422 & vout)
423 
424 !write(std_out,*) 'delocint 08'
425 !##########################################################
426 !### 08. Initialize or update the hessian matrix
427 
428 !Initialise the Hessian matrix using gmet
429  if (itime==1)then
430 
431 !  Initialise the Hessian matrix with ab_mover%userrc.
432 !  this has become unusable because it imposes ndim >= 3 natom
433 !  ident = 3x3 identity matrix
434 !  call hessinit(ab_mover, hessin, gmet, ndim, ucvol)
435    hessin = zero
436    do ii=1, ndim
437      hessin (ii,ii) = one
438    end do
439 
440 !  ! Initialize inverse hessian with identity matrix
441 !  ! in cartesian coordinates, which makes use of metric tensor gmet
442 !  ! in reduced coordinates.
443 !  hessin(:,:)=zero
444 !  do ii=1,ab_mover%natom
445 !  do kk=1,3
446 !  do jj=1,3
447 !  ! Warning : implemented in reduced coordinates
448 !  if (ab_mover%iatfix(kk,ii)==0 .and.&
449 !  & ab_mover%iatfix(jj,ii)==0 )then
450 !  hessin(kk+3*(ii-1),jj+3*(ii-1))=gmet(kk,jj)
451 !  end if
452 !  end do
453 !  end do
454 !  end do
455 !  if(ab_mover%optcell/=0)then
456 !  ! These values might lead to too large changes in some cases
457 !  diag=ab_mover%strprecon*30.0_dp/ucvol
458 !  if(ab_mover%optcell==1) diag=diag/three
459 !  do ii=3*ab_mover%natom+1,ndim
460 !  hessin(ii,ii)=diag
461 !  end do
462 !  end if
463 
464    if (ab_mover%restartxf/=0) then
465 
466      call xfh_recover_deloc(ab_xfh,ab_mover,acell,acell0,cycl_main,&
467 &     residual,hessin,ndim,rprimd,rprimd0,strten,ucvol,ucvol0,vin,vin_prev,&
468 &     vout,vout_prev,xred,deloc,deloc_int,deloc_force,bt_inv_matrix,gprimd,prim_int,&
469 &     u_matrix)
470 
471    end if
472 
473  end if
474 
475  ABI_DEALLOCATE(prim_int)
476 
477  if(itime>1)then
478 !  Update the hessian matrix, by taking into account the
479 !  current pair (x,f) and the previous one.
480    call hessupdt(hessin,ab_mover%iatfix,ab_mover%natom,ndim,vin,&
481 &   vin_prev,vout,vout_prev)
482 
483  end if
484 
485 !DEBUG (vin,vout and hessin before prediction)
486  if(DEBUG)then
487    write(std_out,*) 'Vectors vin and vout and inverse of Hessian (hessin) [before]'
488    write(std_out,*) 'vin:'
489    do ii=1,ndim,3
490      if (ii+2<=ndim)then
491        write(std_out,*) ii,vin(ii:ii+2)
492      else
493        write(std_out,*) ii,vin(ii:ndim)
494      end if
495    end do
496    write(std_out,*) 'vout:'
497    do ii=1,ndim,3
498      if (ii+2<=ndim)then
499        write(std_out,*) ii,vout(ii:ii+2)
500      else
501        write(std_out,*) ii,vout(ii:ndim)
502      end if
503    end do
504    write(std_out,*) 'Inverse Hessian (hessin): ',ndim,'x',ndim
505    do kk=1,ndim
506      do jj=1,ndim,3
507        if (jj+2<=ndim)then
508          write(std_out,*) jj,hessin(jj:jj+2,kk)
509        else
510          write(std_out,*) jj,hessin(jj:ndim,kk)
511        end if
512      end do
513    end do
514  end if
515 
516 !write(std_out,*) 'delocint 09'
517 !##########################################################
518 !### 09. Compute the next values
519 
520  if(ionmov==10 .or. itime==1)then
521 
522 !  Previous cartesian coordinates
523    vin_prev(:)=vin(:)
524 
525 !  New atomic cartesian coordinates are obtained from vin, hessin
526 !  and vout
527    do ii=1,ndim
528      vin(:)=vin(:)-hessin(:,ii)*vout(ii)
529    end do
530 !  Previous atomic forces
531    vout_prev(:)=vout(:)
532 
533  else
534    if(ionmov==11)then
535      ihist_prev = abihist_findIndex(hist,-1)
536      etotal_prev=hist%etot(ihist_prev)
537 !    Here the BFGS algorithm, modified to take into account the
538 !    energy
539      call brdene(etotal,etotal_prev,hessin,&
540 &     ndim,vin,vin_prev,vout,vout_prev)
541 
542    end if
543 
544 !  DEBUG (vin,vout and hessin after prediction)
545    if(DEBUG)then
546      write(std_out,*) 'Vectors vin and vout [after prediction]'
547      write(std_out,*) 'vin:'
548      do ii=1,ndim,3
549        if (ii+2<=ndim)then
550          write(std_out,*) ii,vin(ii:ii+2)
551        else
552          write(std_out,*) ii,vin(ii:ndim)
553        end if
554      end do
555      write(std_out,*) 'vout:'
556      do ii=1,ndim,3
557        if (ii+2<=ndim)then
558          write(std_out,*) ii,vout(ii:ii+2)
559        else
560          write(std_out,*) ii,vout(ii:ndim)
561        end if
562      end do
563    end if
564 
565 !  Implement fixing of atoms : put back old values for fixed
566 !  components
567    do kk=1,ab_mover%natom
568      do jj=1,3
569 !      Warning : implemented in reduced coordinates
570        if ( ab_mover%iatfix(jj,kk)==1) then
571          vin(jj+(kk-1)*3)=vin_prev(jj+(kk-1)*3)
572        end if
573      end do
574    end do
575  end if
576 
577 !write(std_out,*) 'delocint 10'
578 !##########################################################
579 !### 10. Convert from delocalized to xcart and xred
580 
581 !Transfer vin  to deloc_int, acell and rprimd
582  call xfpack_vin2x(acell, acell0, ab_mover%natom-1, ndim,&
583 & ab_mover%nsym, ab_mover%optcell, rprimd, rprimd0,&
584 & ab_mover%symrel, ucvol, ucvol0,&
585 & vin, deloc_int)
586 
587  if(DEBUG)then
588    write (std_out,*) 'Delocalized internals (deloc_int) [after prediction]:'
589    write(std_out,*) 'deloc_int:'
590    do ii=1,ndeloc,3
591      if (ii+2<=ndeloc)then
592        write(std_out,*) ii,deloc_int(ii:ii+2)
593      else
594        write(std_out,*) ii,deloc_int(ii:ndeloc)
595      end if
596    end do
597    write(std_out,*) 'BT Inverse Matrix:'
598    do ii=1,3*ab_mover%natom
599      write(std_out,*) bt_inv_matrix(:,ii)
600    end do
601    write (std_out,*) 'xcart (before deloc2xcart):'
602    do ii=1,ab_mover%natom
603      write (std_out,*) xcart(:,ii)
604    end do
605  end if
606 
607 !this routine contains an iterative scheme to find xcart
608 !from the non-linear relations between deloc and xcart
609 !SIGNIFICANTLY DIFFERENT FROM xcart2deloc
610  call deloc2xcart(deloc,ab_mover%natom,rprimd,xcart,&
611 & deloc_int,bt_inv_matrix,u_matrix)
612 
613 !Convert new xcart (cartesian) to xred (reduced coordinates)
614  call xcart2xred(ab_mover%natom,rprimd,xcart,xred)
615 
616 
617 !write(std_out,*) 'delocint 11'
618 !##########################################################
619 !### 11. Update the history with the prediction
620 
621 !Increase indexes
622  hist%ihist = abihist_findIndex(hist,+1)
623 
624  if(ab_mover%optcell/=0)then
625    call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
626  end if
627 
628 !Fill the history with the variables
629 !xred, acell, rprimd, vel
630  call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
631  ihist_prev = abihist_findIndex(hist,-1)
632  hist%vel(:,:,hist%ihist)=hist%vel(:,:,ihist_prev)
633 
634  if(zDEBUG)then
635    write (std_out,*) 'residual:'
636    do kk=1,ab_mover%natom
637      write (std_out,*) residual(:,kk)
638    end do
639    write (std_out,*) 'strten:'
640    write (std_out,*) strten(1:3),ch10,strten(4:6)
641    write (std_out,*) 'etotal:'
642    write (std_out,*) etotal
643  end if
644 
645 end subroutine pred_delocint

ABINIT/xcart2deloc [ Functions ]

[ Top ] [ 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

PARENTS

      deloc2xcart,pred_delocint,xfh_recover_deloc

CHILDREN

SOURCE

1585 subroutine xcart2deloc(deloc,natom,rprimd,xcart,bt_inv_matrix,u_matrix,deloc_int,prim_int)
1586 
1587 
1588 !This section has been created automatically by the script Abilint (TD).
1589 !Do not modify the following lines by hand.
1590 #undef ABI_FUNC
1591 #define ABI_FUNC 'xcart2deloc'
1592 !End of the abilint section
1593 
1594  implicit none
1595 
1596 !Arguments ------------------------------------
1597 !scalars
1598  integer,intent(in) :: natom
1599  type(delocint),intent(in) :: deloc
1600 !arrays
1601  real(dp),intent(in) :: rprimd(3,3),xcart(3,natom)
1602  real(dp),intent(inout) :: u_matrix(deloc%ninternal,3*(natom-1))
1603  real(dp),intent(out) :: bt_inv_matrix(3*(natom-1),3*natom)
1604  real(dp),intent(out) :: deloc_int(3*(natom-1))
1605  real(dp),intent(out) :: prim_int(deloc%ninternal)
1606 
1607 !Local variables-------------------------------
1608 !scalars
1609 integer :: ii
1610 logical :: DEBUG=.FALSE.
1611 !arrays
1612  real(dp) :: b_matrix(deloc%ninternal,3*natom)
1613 
1614 ! ******************************************************************
1615 
1616  call calc_prim_int(deloc,natom,rprimd,xcart,prim_int)
1617  if (DEBUG)then
1618    write(std_out,*) 'Primitive Internals'
1619    do ii=1,deloc%ninternal
1620      write(std_out,*) prim_int(ii)
1621    end do
1622  end if
1623 
1624  call calc_b_matrix(deloc,natom,rprimd,xcart,b_matrix)
1625  if (DEBUG)then
1626    write(std_out,*) 'B Matrix'
1627    do ii=1,deloc%ninternal
1628      write(std_out,*) b_matrix(:,ii)
1629    end do
1630  end if
1631 
1632  call calc_btinv_matrix(b_matrix,natom,deloc%ninternal,&
1633 & bt_inv_matrix,u_matrix)
1634  if (DEBUG)then
1635    write(std_out,*) 'BT Inverse Matrix'
1636    do ii=1,3*natom
1637      write(std_out,*) bt_inv_matrix(:,ii)
1638    end do
1639  end if
1640 
1641 !calculate value of delocalized internals
1642 
1643  call dgemv('T',deloc%ninternal,3*(natom-1),one,&
1644 & u_matrix,deloc%ninternal,prim_int,1,zero,deloc_int,1)
1645 
1646 end subroutine xcart2deloc

ABINIT/xfh_recover_deloc [ Functions ]

[ Top ] [ Functions ]

NAME

 xfh_recover_deloc

FUNCTION

 Update the contents of the history xfhist taking values
 from xred, acell, rprim, fred_corrected and strten

INPUTS

OUTPUT

PARENTS

      pred_delocint

CHILDREN

      fred2fdeloc,hessupdt,xcart2deloc,xfpack_f2vout,xfpack_x2vin,xred2xcart

SOURCE

1868 subroutine xfh_recover_deloc(ab_xfh,ab_mover,acell,acell0,cycl_main,&
1869 & fred,hessin,ndim,rprim,rprimd0,strten,ucvol,ucvol0,vin,vin_prev,&
1870 & vout,vout_prev,xred,deloc,deloc_int,deloc_force,btinv,gprimd,prim_int,&
1871 & u_matrix)
1872 
1873 
1874 !This section has been created automatically by the script Abilint (TD).
1875 !Do not modify the following lines by hand.
1876 #undef ABI_FUNC
1877 #define ABI_FUNC 'xfh_recover_deloc'
1878 !End of the abilint section
1879 
1880 implicit none
1881 
1882 !Arguments ------------------------------------
1883 !scalars
1884 
1885 integer,intent(in) :: ndim
1886 integer,intent(out) :: cycl_main
1887 real(dp),intent(inout) :: ucvol,ucvol0
1888 type(ab_xfh_type),intent(inout) :: ab_xfh
1889 type(abimover),intent(in) :: ab_mover
1890 ! DELOCINT specials
1891 type(delocint),intent(in) :: deloc
1892 
1893 !arrays
1894 real(dp),intent(inout) :: acell(3)
1895 real(dp),intent(in) :: acell0(3)
1896 real(dp),intent(inout) :: hessin(:,:)
1897 real(dp),intent(inout) :: xred(3,ab_mover%natom)
1898 real(dp),intent(inout) :: rprim(3,3)
1899 real(dp),intent(inout) :: rprimd0(3,3)
1900 real(dp),intent(inout) :: fred(3,ab_mover%natom)
1901 real(dp),intent(inout) :: strten(6)
1902 real(dp),intent(inout) :: vin(:)
1903 real(dp),intent(inout) :: vin_prev(:)
1904 real(dp),intent(inout) :: vout(:)
1905 real(dp),intent(inout) :: vout_prev(:)
1906 ! DELOCINT specials
1907 real(dp),intent(inout) :: deloc_force(3*(ab_mover%natom-1))
1908 real(dp),intent(inout) :: deloc_int(3*(ab_mover%natom-1))
1909 real(dp),intent(inout) :: btinv(3*(ab_mover%natom-1),3*ab_mover%natom)
1910 real(dp),intent(inout) :: prim_int(:),u_matrix(:,:),gprimd(3,3)
1911 
1912 !Local variables-------------------------------
1913 !scalars
1914 integer :: ixfh
1915 real(dp) :: xcart(3,ab_mover%natom)
1916 
1917 !*********************************************************************
1918 
1919  if(ab_xfh%nxfh/=0)then
1920 !  Loop over previous time steps
1921    do ixfh=1,ab_xfh%nxfh
1922 
1923 !    For that time step, get new (x,f) from xfhist
1924      xred(:,:)     =ab_xfh%xfhist(:,1:ab_mover%natom        ,1,ixfh)
1925      rprim(1:3,1:3)=ab_xfh%xfhist(:,ab_mover%natom+2:ab_mover%natom+4,1,ixfh)
1926      acell(:)      =ab_xfh%xfhist(:,ab_mover%natom+1,1,ixfh)
1927      fred(:,:)     =ab_xfh%xfhist(:,1:ab_mover%natom,2,ixfh)
1928 !    This use of results_gs is unusual
1929      strten(1:3)   =ab_xfh%xfhist(:,ab_mover%natom+2,2,ixfh)
1930      strten(4:6)   =ab_xfh%xfhist(:,ab_mover%natom+3,2,ixfh)
1931 
1932 !    !DEBUG
1933 !    write (ab_out,*) '---READ FROM XFHIST---'
1934 
1935 !    write (ab_out,*) 'XRED'
1936 !    do kk=1,ab_mover%natom
1937 !    write (ab_out,*) xred(:,kk)
1938 !    end do
1939 !    write (ab_out,*) 'FRED'
1940 !    do kk=1,ab_mover%natom
1941 !    write (ab_out,*) fred(:,kk)
1942 !    end do
1943 !    write(ab_out,*) 'RPRIM'
1944 !    do kk=1,3
1945 !    write(ab_out,*) rprim(:,kk)
1946 !    end do
1947 !    write(ab_out,*) 'ACELL'
1948 !    write(ab_out,*) acell(:)
1949 !    !DEBUG
1950 
1951 !    Convert input xred (reduced coordinates) to xcart (cartesian)
1952      call xred2xcart(ab_mover%natom,rprimd0,xcart,xred)
1953 !    Convert input coordinates in Delocalized internals
1954      call xcart2deloc(deloc,ab_mover%natom,rprimd0,xcart,&
1955 &     btinv,u_matrix,deloc_int,prim_int)
1956 !    Convert forces to delocalized coordinates for next step
1957      call fred2fdeloc(btinv,deloc_force,fred,ab_mover%natom,gprimd)
1958 
1959 !    Transfer it in vin, vout
1960      call xfpack_x2vin(acell,acell0,ab_mover%natom-1,&
1961 &     ndim,ab_mover%nsym,ab_mover%optcell,rprim,rprimd0,&
1962 &     ab_mover%symrel,ucvol,ucvol0,vin,deloc_int)
1963      call xfpack_f2vout(deloc_force,ab_mover%natom-1,&
1964 &     ndim,ab_mover%optcell,ab_mover%strtarget,strten,&
1965 &     ucvol,vout)
1966 !    Get old time step, if any, and update inverse hessian
1967      if(ixfh/=1)then
1968        xred(:,:)     =ab_xfh%xfhist(:,1:ab_mover%natom,1,ixfh-1)
1969        rprim(1:3,1:3)=&
1970 &       ab_xfh%xfhist(:,ab_mover%natom+2:ab_mover%natom+4,1,ixfh-1)
1971        acell(:)=ab_xfh%xfhist(:,ab_mover%natom+1,1,ixfh-1)
1972        fred(:,:)=ab_xfh%xfhist(:,1:ab_mover%natom,2,ixfh-1)
1973 !      This use of results_gs is unusual
1974        strten(1:3)=ab_xfh%xfhist(:,ab_mover%natom+2,2,ixfh-1)
1975        strten(4:6)=ab_xfh%xfhist(:,ab_mover%natom+3,2,ixfh-1)
1976 
1977 !      Convert input xred (reduced coordinates) to xcart (cartesian)
1978        call xred2xcart(ab_mover%natom,rprimd0,xcart,xred)
1979 !      Convert input coordinates in Delocalized internals
1980        call xcart2deloc(deloc,ab_mover%natom,rprimd0,xcart,&
1981 &       btinv,u_matrix,deloc_int,prim_int)
1982 !      Convert forces to delocalized coordinates for next step
1983        call fred2fdeloc(btinv,deloc_force,fred,ab_mover%natom,gprimd)
1984 
1985 !      Tranfer it in vin_prev, vout_prev
1986        call xfpack_x2vin(acell,acell0,ab_mover%natom-1,&
1987 &       ndim,ab_mover%nsym,ab_mover%optcell,rprim,rprimd0,&
1988 &       ab_mover%symrel,ucvol,ucvol0,vin_prev,deloc_int)
1989        call xfpack_f2vout(deloc_force,ab_mover%natom-1,&
1990 &       ndim,ab_mover%optcell,ab_mover%strtarget,strten,&
1991 &       ucvol,vout_prev)
1992 
1993 !      write(ab_out,*) 'Hessian matrix before update',ndim,'x',ndim
1994 !      write(ab_out,*) 'ixfh=',ixfh
1995 !      do kk=1,ndim
1996 !      do jj=1,ndim,3
1997 !      if (jj+2<=ndim)then
1998 !      write(ab_out,*) jj,hessin(jj:jj+2,kk)
1999 !      else
2000 !      write(ab_out,*) jj,hessin(jj:ndim,kk)
2001 !      end if
2002 !      end do
2003 !      end do
2004 
2005        call hessupdt(hessin,ab_mover%iatfix,ab_mover%natom-1,ndim,&
2006 &       vin,vin_prev,vout,vout_prev)
2007 
2008 !      !DEBUG
2009 !      write(ab_out,*) 'Hessian matrix after update',ndim,'x',ndim
2010 !      do kk=1,ndim
2011 !      do jj=1,ndim,3
2012 !      if (jj+2<=ndim)then
2013 !      write(ab_out,*) jj,hessin(jj:jj+2,kk)
2014 !      else
2015 !      write(ab_out,*) jj,hessin(jj:ndim,kk)
2016 !      end if
2017 !      end do
2018 !      end do
2019 !      !DEBUG
2020 
2021      end if !if(ab_xfh%nxfh/=0)
2022 
2023 !    End loop over previous time steps
2024    end do
2025 
2026 !  The hessian has been generated,
2027 !  as well as the latest vin and vout
2028 !  so will cycle the main loop
2029    cycl_main=1
2030 
2031  end if
2032 
2033 end subroutine xfh_recover_deloc