Révision 4
| utils/AnaPathrel (revision 4) | ||
|---|---|---|
| 1 | 1 |
#!/bin/bash |
| 2 |
if [ $# -lt 4 ]; then
|
|
| 3 |
echo "Use: $0 File.out MaxCyc PathName NGeomF"
|
|
| 2 |
if [ $# -lt 1 ]; then
|
|
| 3 |
echo "Use: $0 File.out" |
|
| 4 | 4 |
exit |
| 5 | 5 |
fi |
| 6 | 6 |
|
| 7 | 7 |
Fout=$1 |
| 8 |
ItMax=$2
|
|
| 9 |
Nom=$3
|
|
| 10 |
NGeomF=$4
|
|
| 8 |
ItMax=`grep MAXCYC $Fout | tail -1 | awk '{print $3}'`
|
|
| 9 |
Nom=`grep PATHNAME $Fout | tail -1 | awk '{print $3}'`
|
|
| 10 |
NGeomF=`grep NGEOMF $Fout | tail -1 | awk '{print $3}'`
|
|
| 11 | 11 |
|
| 12 |
|
|
| 12 | 13 |
export LANG=C |
| 13 | 14 |
|
| 14 | 15 |
echo "#ItMax=$ItMax" |
| utils/Xyz2Path.param (revision 4) | ||
|---|---|---|
| 1 |
INTEGER*4 MaxNAt |
|
| 2 |
PARAMETER (MaxNat=500) |
|
| 1 |
INTEGER(4), PARAMETER :: MaxNAt=500 |
|
| 2 |
|
|
| utils/AnaPath (revision 4) | ||
|---|---|---|
| 1 | 1 |
#!/bin/bash |
| 2 |
if [ $# -lt 4 ]; then
|
|
| 3 |
echo "Use: $0 File.out MaxCyc PathName NGeomF [x column]"
|
|
| 2 |
if [ $# -lt 1 ]; then
|
|
| 3 |
echo "Use: $0 File.out [x column]" |
|
| 4 | 4 |
exit |
| 5 | 5 |
fi |
| 6 | 6 |
|
| 7 | 7 |
Fout=$1 |
| 8 |
ItMax=$2 |
|
| 9 |
Nom=$3 |
|
| 10 |
NGeomF=$4 |
|
| 11 | 8 |
|
| 9 |
ItMax=`grep MAXCYC $Fout | tail -1 | awk '{print $3}'`
|
|
| 10 |
Nom=`grep PATHNAME $Fout | tail -1 | awk '{print $3}'`
|
|
| 11 |
NGeomF=`grep NGEOMF $Fout | tail -1 | awk '{print $3}'`
|
|
| 12 |
|
|
| 12 | 13 |
xcol=1 |
| 13 |
if [ $# -eq 5 ]; then
|
|
| 14 |
xcol=$5
|
|
| 14 |
if [ $# -eq 2 ]; then
|
|
| 15 |
xcol=$2
|
|
| 15 | 16 |
fi |
| 16 | 17 |
|
| 17 | 18 |
export LANG=C |
| ... | ... | |
| 41 | 42 |
fi |
| 42 | 43 |
echo "Creating $Nom.datl" |
| 43 | 44 |
|
| 45 |
ItDone=-1 |
|
| 44 | 46 |
for i in `seq 0 $ItMax` |
| 45 | 47 |
do |
| 46 |
xyz2path ${Nom}${Ext}.$i
|
|
| 47 |
cat Scan.dat >> $Nom.datl |
|
| 48 |
echo " " >> $Nom.datl |
|
| 49 |
echo " " >> $Nom.datl |
|
| 48 |
if [ -s ${Nom}${Ext}.$i ]; then
|
|
| 49 |
xyz2path ${Nom}${Ext}.$i
|
|
| 50 |
cat Scan.dat >> $Nom.datl |
|
| 51 |
echo " " >> $Nom.datl |
|
| 52 |
echo " " >> $Nom.datl |
|
| 53 |
let ItDone=ItDone+1 |
|
| 54 |
fi |
|
| 50 | 55 |
done |
| 51 | 56 |
|
| 52 | 57 |
icol=`tail -1 Scan.dat | wc -w` |
| ... | ... | |
| 58 | 63 |
EOF |
| 59 | 64 |
|
| 60 | 65 |
|
| 61 |
for i in `seq 1 $ItMax`
|
|
| 66 |
for i in `seq 1 $ItDone`
|
|
| 62 | 67 |
do |
| 63 | 68 |
echo "plot \"$Nom.datl\" i 0 u xcol:$icol w lp " >> ${Nom}_l.gplot
|
| 64 | 69 |
echo "replot \"$Nom.datl\" i $i u xcol:$icol w lp " >> ${Nom}_l.gplot
|
| ... | ... | |
| 76 | 81 |
EOF |
| 77 | 82 |
|
| 78 | 83 |
echo "plot \"$Nom.datl\" i 0 u xcol:$icol w lp " >> ${Nom}_l2.gplot
|
| 79 |
for i in `seq 1 $ItMax`
|
|
| 84 |
for i in `seq 1 $ItDone`
|
|
| 80 | 85 |
do |
| 81 | 86 |
echo "replot \"$Nom.datl\" i $i u xcol:$icol w lp " >> ${Nom}_l2.gplot
|
| 82 | 87 |
done |
| ... | ... | |
| 89 | 94 |
EOF |
| 90 | 95 |
|
| 91 | 96 |
echo "plot \"$Nom.datl\" i 0 u 0:$icol w lp " >> ${Nom}_l3.gplot
|
| 92 |
for i in `seq 1 $ItMax`
|
|
| 97 |
for i in `seq 1 $ItDone`
|
|
| 93 | 98 |
do |
| 94 | 99 |
echo "replot \"$Nom.datl\" i $i u 0:$icol w lp " >> ${Nom}_l3.gplot
|
| 95 | 100 |
done |
| doc/Mini_help.tex (revision 4) | ||
|---|---|---|
| 191 | 191 |
\item Step\_method: method to compute the step for optimizing a geometry; choices are: |
| 192 | 192 |
- RFO: Rational function optimization |
| 193 | 193 |
- GDIIS: Geometry optimization using direct inversion in the iterative supspace |
| 194 |
\item HesUpd: method to update the hessian. By default, it is Murtagh-Sargent
|
|
| 194 |
\item HUpdate: method to update the hessian. By default, it is Murtagh-Sargent
|
|
| 195 | 195 |
Except for geometry optimization where it is BFGS. |
| 196 | 196 |
\item MaxCyc: maximum number of iterations for the path optimization |
| 197 | 197 |
\item Smax: Maximum length of a step during path optimization |
| ... | ... | |
| 217 | 217 |
Flags: |
| 218 | 218 |
\begin{description}
|
| 219 | 219 |
\item MW: Flag. True if one wants to work in Mass Weighted coordinates. Default=.TRUE. |
| 220 |
\item Renum: Flag. True if one wants to reoder the atoms in the initial order. default is .TRUE. most of the time. |
|
| 220 |
\item Renum: Flag. True if one wants to reoder the atoms in the |
|
| 221 |
initial order. default is .TRUE. most of the time. |
|
| 221 | 222 |
\item OptProd: True if one wants to optimize the geometry of the products. |
| 222 | 223 |
\item OptReac: True if one wants to optimize the geometry of the reactants. |
| 224 |
\item CalcEProd: if TRUE the product energy will be |
|
| 225 |
computed. Default is FALSE. Not Compatible with RunMode=Para". |
|
| 226 |
\item CalcEReac: if TRUE the reactants energy will be |
|
| 227 |
computed. Default is FALSE. Not Compatible with RunMode=Para". |
|
| 223 | 228 |
\item PathOnly:TRUE if one wants to generate the initial path, and stops. |
| 224 | 229 |
\item Hinv: if True, then Hessian inversed is used. |
| 225 |
\item IniHup: if True, then Hessian inverse is extrapolated using the initial path calculations. |
|
| 226 |
\item HupNeighbour: if True, then Hessian inverse is extrapolated using the neighbouring points of the path. |
|
| 230 |
\item IniHup: if True, then Hessian inverse is extrapolated |
|
| 231 |
using the initial path calculations. |
|
| 232 |
\item HupNeighbour: if True, then Hessian inverse is |
|
| 233 |
extrapolated using the neighbouring points of the path. |
|
| 227 | 234 |
\item FFrozen: True if one wants to freeze the positions of some atoms. |
| 228 |
If True, a \verb!&frozenlist! namelist containing the list of frozen atoms must be given. |
|
| 235 |
If True, a \verb!&frozenlist! namelist containing the |
|
| 236 |
list of frozen atoms must be given. |
|
| 229 | 237 |
If VASP is used, and frozen is not given |
| 230 | 238 |
here, the program will use the F flags of the input geometry |
| 231 | 239 |
\item FCart: True if one wants to describe some atoms using cartesian coordinates. |
| ... | ... | |
| 251 | 259 |
Serial or Parallel. |
| 252 | 260 |
\begin{itemize}
|
| 253 | 261 |
\item[Serial] |
| 254 |
In the Serial mode, Path uses Vasp to compute the energy and forces of each image separately (as it does for Gaussian for example). |
|
| 255 |
Therefore, the INCAR file should not contain any references to the number of images (no IMAGES command!). |
|
| 262 |
In the Serial mode, Path uses Vasp to compute the energy and forces of |
|
| 263 |
each image separately (as it does for Gaussian for example). |
|
| 264 |
Therefore, the INCAR file should not contain any references to the |
|
| 265 |
number of images (no IMAGES command!). |
|
| 256 | 266 |
|
| 257 | 267 |
\item[Parallel] |
| 258 |
In the Parallel mode, Path uses VASP to compute at once the energies and forces for all the images (as VASP does for a NEB calculation). |
|
| 268 |
In the Parallel mode, Path uses VASP to compute at once the energies |
|
| 269 |
and forces for all the images (as VASP does for a NEB calculation). |
|
| 259 | 270 |
Therefore, the INCAR file MUST contain the IMAGES command. |
| 260 | 271 |
More, the directories 00, 01, ... should exist. |
| 261 | 272 |
|
| src/Egrad.f90 (revision 4) | ||
|---|---|---|
| 6 | 6 |
|
| 7 | 7 |
use Path_module, only : Nat,AtName,Coord,dzdc,indzmat,Nom,Atome,massat,unit, & |
| 8 | 8 |
prog,NCart,XyzGeomF,IntCoordF,XyzGeomI, renum, OrderInv |
| 9 |
! IntCoordF(NGeomF,NCoord),GeomOld_all(NGeomF,NCoord)
|
|
| 10 |
! allocated in Path.f90
|
|
| 9 |
! IntCoordF(NGeomF,NCoord),GeomOld_all(NGeomF,NCoord)
|
|
| 10 |
! allocated in Path.f90
|
|
| 11 | 11 |
use Io_module |
| 12 | 12 |
IMPLICIT NONE |
| 13 | 13 |
|
| ... | ... | |
| 161 | 161 |
Idx=Idx+3 |
| 162 | 162 |
END DO |
| 163 | 163 |
DEALLOCATE(GradTmp) |
| 164 |
CASE ('BAKER')
|
|
| 164 |
CASE ('BAKER')
|
|
| 165 | 165 |
WRITE (*,*) "Coord=",TRIM(COORD),": use Egrad_baker.f90, instead of Egrad.f90. Stop" |
| 166 | 166 |
CASE ('MIXED')
|
| 167 | 167 |
ALLOCATE(GradTmp(3*Nat)) |
| src/Hupdate_all.f90 (revision 4) | ||
|---|---|---|
| 1 | 1 |
subroutine hupdate_all(n, ds, dgrad, Hess) |
| 2 | 2 |
|
| 3 |
use Path_module, only : hesupd,Hinv
|
|
| 3 |
use Path_module, only : HUpdate,Hinv
|
|
| 4 | 4 |
use Io_module, only : IoOut |
| 5 | 5 |
|
| 6 | 6 |
IMPLICIT NONE |
| ... | ... | |
| 14 | 14 |
|
| 15 | 15 |
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
| 16 | 16 |
! |
| 17 |
! This subroutine is a drive for the update of the Hessian |
|
| 17 |
! This subroutine is a driver for the update of the Hessian
|
|
| 18 | 18 |
! in the Quasi-Newton optimization step. |
| 19 | 19 |
! We propose mostly update of the inverse Hessian. |
| 20 | 20 |
! |
| ... | ... | |
| 29 | 29 |
!!!!!!!!!!!!!!!! |
| 30 | 30 |
! |
| 31 | 31 |
! For now, we propose : |
| 32 |
! - For update of inverser Hessian: |
|
| 33 |
! - BFGS |
|
| 32 |
! - For update of inverse Hessian: |
|
| 33 |
! - BFGS (Broyden, Fletcher, Goldfarb, and Shanno) |
|
| 34 |
! - DFP (Davidson - Fletcher - Powell ) |
|
| 35 |
! - MS (Murtagh-Sargent) |
|
| 34 | 36 |
! |
| 35 | 37 |
! - For update of hessian: |
| 36 | 38 |
! - Bofill (combination of Murtagh-Sargent and Powell-symmetric-Broyden |
| 37 |
! - Murtagh-Sargent |
|
| 39 |
! - MS (Murtagh-Sargent) |
|
| 40 |
! - BFGS |
|
| 41 |
! - DFP |
|
| 42 |
!---------------------------------------------------- |
|
| 43 |
! We use the fact that some formula are the dual of other |
|
| 44 |
! that is, one can convert the formula for update of the Hessian |
|
| 45 |
! into the update of the inverse Hessian by replacing the Hessian |
|
| 46 |
! (denoted by B in the Fletcher book and in the Nocedal&Wrigth book) |
|
| 47 |
! by the inverse hessian (denoted by H), and replacing ds by dgrad in the |
|
| 48 |
! formula. |
|
| 49 |
! See for example: |
|
| 50 |
! 1) Numerical Optimization, Springer 2nd Ed., 2006 |
|
| 51 |
! by Jorge Nocedal & Stephen J. Wright |
|
| 52 |
! 2) ractical Methods of Optimization, John Wiley & Sons, second ed., 1987. |
|
| 53 |
! by R. Fletcher |
|
| 38 | 54 |
! |
| 55 |
! |
|
| 39 | 56 |
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
| 40 | 57 |
|
| 41 | 58 |
|
| ... | ... | |
| 53 | 70 |
|
| 54 | 71 |
! ====================================================================== |
| 55 | 72 |
|
| 56 |
debug = valid ('DEBUG HESS') .OR. Valid('DEBUG HUPDATE')
|
|
| 73 |
debug = valid ('HESS') .OR. Valid('HUPDATE')
|
|
| 57 | 74 |
|
| 58 | 75 |
if (debug) WRITE(*,*) "================================= Entering Hupdate_all ===============" |
| 59 | 76 |
|
| 60 | 77 |
SELECT CASE (Hinv) |
| 61 | 78 |
CASE (.TRUE.) |
| 62 |
SELECT CASE (HesUpd)
|
|
| 79 |
SELECT CASE (HUpdate)
|
|
| 63 | 80 |
CASE ("BFGS")
|
| 64 | 81 |
Call Hinvup_BFGS(n,ds,dgrad,hess) |
| 82 |
CASE ("MS")
|
|
| 83 |
! we use the fact that MS is self-dual |
|
| 84 |
Call Hupdate_MS(n,dgrad,ds,hess) |
|
| 85 |
CASE ("DFP")
|
|
| 86 |
Call Hinvup_DFP(n,ds,dgrad,hess) |
|
| 65 | 87 |
CASE DEFAULT |
| 66 |
WRITE(*,*) Trim(HesUpd) // " not known: using BFGS"
|
|
| 88 |
WRITE(*,*) Trim(HUpdate) // " not known: using BFGS"
|
|
| 67 | 89 |
Call Hinvup_BFGS(n,ds,dgrad,hess) |
| 68 | 90 |
END SELECT |
| 69 | 91 |
CASE (.FALSE.) |
| 70 |
SELECT CASE (HesUpd)
|
|
| 92 |
SELECT CASE (HUpdate)
|
|
| 71 | 93 |
CASE ("BOFILL")
|
| 72 | 94 |
Call Hupdate_Bofill(n,ds,dgrad,hess) |
| 73 | 95 |
CASE ("MS")
|
| 74 | 96 |
Call Hupdate_MS(n,ds,dgrad,hess) |
| 97 |
CASE ("DFP")
|
|
| 98 |
! we use the fact that DFP is the dual of BFGS |
|
| 99 |
Call Hinvup_BFGS(n,dgrad,ds,hess) |
|
| 100 |
CASE ("BFGS")
|
|
| 101 |
! we use the fact that DFP is the dual of BFGS |
|
| 102 |
Call Hinvup_DFP(n,dgrad,ds,hess) |
|
| 75 | 103 |
CASE DEFAULT |
| 76 |
WRITE(*,*) Trim(HesUpd) // " not known: using Bofill"
|
|
| 104 |
WRITE(*,*) Trim(HUpdate) // " not known: using Bofill"
|
|
| 77 | 105 |
Call Hupdate_Bofill(n,ds,dgrad,hess) |
| 78 | 106 |
END SELECT |
| 79 | 107 |
END SELECT |
| 80 | 108 |
|
| 109 |
if (debug) WRITE(*,*) "================================= Exiting Hupdate_all ===============" |
|
| 81 | 110 |
end subroutine hupdate_all |
| src/Hupdate_MS.f90 (revision 4) | ||
|---|---|---|
| 1 |
! This program reads an old hessian, old forces, new forces |
|
| 2 |
! and print the updated hessian and old coordinates. |
|
| 3 |
! We use the Murtagh-Sargent update |
|
| 1 |
! Hessian - or its inverse - update |
|
| 2 |
! We use the Murtagh-Sargent update, also called Symmetric Rank 1 |
|
| 3 |
! |
|
| 4 |
! input: |
|
| 5 |
! - nfree: number of coordinates |
|
| 6 |
! - ds: GeomOld-Geom (expressed in the optimization coord.) |
|
| 7 |
! - dgrad: GradOld-Grad (expressed in the optimization coord.) |
|
| 8 |
! |
|
| 9 |
! input/output: |
|
| 10 |
! - Hess: old Hessian in input, replaced by updated Hessian on output |
|
| 11 |
! (or inverse Hessian) |
|
| 12 |
! |
|
| 13 |
! |
|
| 14 |
!!!!!!!!!!!!!!!!!!!!!! |
|
| 15 |
! The update formula for the Hessian is: |
|
| 16 |
! (following Fletcher, approximate Hessian is denoted by B) |
|
| 17 |
! B(k+1)=B+[(y-Bs){(y-Bs)T}]/[{(y-Bs)T}s]
|
|
| 18 |
! |
|
| 19 |
! it is self dual, ie for the inverse Hessian, denote by H: |
|
| 20 |
! H(k+1)=H+[(s-Hy){(s-Hy)T}]/[{(s-Hy)T}y]
|
|
| 21 |
! |
|
| 22 |
! |
|
| 23 |
! We decompose the computation in two steps: |
|
| 24 |
! first we compute the update, and then we add it to Hess to get the |
|
| 25 |
! new Hessian. |
|
| 26 |
! This allows us to introduce a scaling factor if needed later. |
|
| 27 |
! |
|
| 28 |
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|
| 4 | 29 |
|
| 5 | 30 |
SUBROUTINE Hupdate_MS(Nfree,ds,dgrad,Hess) |
| 6 | 31 |
|
| ... | ... | |
| 23 | 48 |
CHARACTER(120) :: fmt,fmt2,Line |
| 24 | 49 |
INTEGER(KINT) :: i,j, idx |
| 25 | 50 |
|
| 26 |
LOGICAL debug,valid |
|
| 27 |
EXTERNAL valid |
|
| 51 |
LOGICAL :: debug |
|
| 28 | 52 |
|
| 53 |
interface |
|
| 54 |
function valid(string) result (isValid) |
|
| 55 |
logical :: isValid |
|
| 56 |
character(*), intent(in) :: string |
|
| 57 |
end function valid |
|
| 58 |
|
|
| 59 |
end interface |
|
| 60 |
|
|
| 61 |
|
|
| 62 |
|
|
| 63 |
|
|
| 29 | 64 |
debug=valid("Hupdate").OR.Valid("Hupdate_MS")
|
| 30 | 65 |
|
| 31 | 66 |
if (debug) Call Header("Entering Hupdate_MS")
|
| 32 | 67 |
|
| 33 |
ALLOCATE(HessOld(nfree,nfree), HupMS(nfree,nfree), &
|
|
| 68 |
ALLOCATE(HessOld(nfree,nfree), & |
|
| 34 | 69 |
Hup(nfree,nfree), & |
| 35 | 70 |
HupInt(nfree), HS(nfree), & |
| 36 | 71 |
Num(nfree,nfree)) |
| ... | ... | |
| 58 | 93 |
HS(I)=HS(I)+HessOld(i,j)*Ds(j) |
| 59 | 94 |
END DO |
| 60 | 95 |
END DO |
| 96 |
|
|
| 97 |
if (debug) WRITE(*,*) "HS:",HS |
|
| 61 | 98 |
HupInt=DGrad-HS |
| 62 | 99 |
|
| 63 | 100 |
! WE calculate the MS update |
| ... | ... | |
| 69 | 106 |
END DO |
| 70 | 107 |
Den=Den+Hupint(i)*Ds(i) |
| 71 | 108 |
END DO |
| 109 |
if (debug) WRITE(*,*) "Den:",Den,1./Den |
|
| 72 | 110 |
|
| 73 |
HupMS=Num/Den
|
|
| 111 |
Hup=Num/Den |
|
| 74 | 112 |
|
| 75 |
Hess=HessOld+HupMS
|
|
| 113 |
Hess=HessOld+Hup |
|
| 76 | 114 |
ELSE |
| 77 | 115 |
Hess=HessOld |
| 78 | 116 |
END IF |
| ... | ... | |
| 81 | 119 |
! WRITE(*,fmt) Hess(i,:) |
| 82 | 120 |
! END DO |
| 83 | 121 |
|
| 84 |
DEALLOCATE(HessOld, HupMS,HupInt, HS, Num )
|
|
| 122 |
DEALLOCATE(HessOld,Hup, HupInt, HS, Num )
|
|
| 85 | 123 |
|
| 86 | 124 |
if (debug) Call Header("Hupdate_MS Over")
|
| 87 | 125 |
|
| src/Path_module.f90 (revision 4) | ||
|---|---|---|
| 56 | 56 |
|
| 57 | 57 |
|
| 58 | 58 |
LOGICAL :: Freq, MW, Bohr, Renum, Hinv |
| 59 |
LOGICAL :: OptReac, OptProd, PathOnly |
|
| 59 |
! OptReac: if TRUE the reactants will be optimized. Default is FALSE |
|
| 60 |
LOGICAL :: OptReac |
|
| 61 |
! OptProd: if TRUE the products will be optimized. Default is FALSE |
|
| 62 |
LOGICAL :: OptProd |
|
| 63 |
! PathOnly: if TRUE, only the initial path will be constructed. |
|
| 64 |
! Default is FALSE, that is OpenPath construct the path and optimize it |
|
| 65 |
LOGICAL :: PathOnly |
|
| 66 |
! CalcEReac: if TRUE the reactants energy will be computed. Default is FALSE |
|
| 67 |
! Not compatible with RunMode=Para |
|
| 68 |
LOGICAL :: CalcEReac |
|
| 69 |
! CalEProd: if TRUE the products energy will be computed. Default is FALSE |
|
| 70 |
! Not compatible with RunMode=Para |
|
| 71 |
LOGICAL :: CalcEProd |
|
| 60 | 72 |
LOGICAL :: IniHup,HupNeighbour |
| 61 | 73 |
|
| 62 | 74 |
REAL(KREAL) :: Fact,FTan |
| ... | ... | |
| 79 | 91 |
CHARACTER(SCHARS) :: Prog="GAUSSIAN" |
| 80 | 92 |
CHARACTER(VLCHARS) :: ProgExe="g03" |
| 81 | 93 |
|
| 82 |
CHARACTER(SCHARS) :: runtyp='GEOMETRY OPTIMIZATION' |
|
| 94 |
! PFL 06/2011: HesUpd is deprecated |
|
| 95 |
! I replace it by HUpdate that is more intuitive to me |
|
| 83 | 96 |
CHARACTER(SCHARS) :: hesupd="BFGS" |
| 97 |
! HUpdate indicates which method to use for updating the Hessian or its inverse |
|
| 98 |
CHARACTER(SCHARS) :: HUpdate="BFGS" |
|
| 84 | 99 |
|
| 85 | 100 |
REAL(KREAL), ALLOCATABLE, TARGET :: XyzGeomI(:,:,:) ! (NGeomI,3,Nat) |
| 86 | 101 |
REAL(KREAL), ALLOCATABLE :: XyzGeomF(:,:,:) ! (NGeomF,3,Nat) |
| src/Extrapol_baker.f90 (revision 4) | ||
|---|---|---|
| 166 | 166 |
|
| 167 | 167 |
|
| 168 | 168 |
call CalcRmsd(Nat,XyzTmp(1:Nat,1),XyzTmp(1:Nat,2),XyzTmp(1:Nat,3), & |
| 169 |
XyzTmp2(1:Nat,1),XyzTmp2(1:Nat,2),XyzTmp2(1:Nat,3), &
|
|
| 169 |
XyzTmp2(1:Nat,1),XyzTmp2(1:Nat,2),XyzTmp2(1:Nat,3), &
|
|
| 170 | 170 |
MRot,rmsd,.TRUE.,.TRUE.) |
| 171 | 171 |
|
| 172 | 172 |
|
| src/Extrapol_int.f90 (revision 4) | ||
|---|---|---|
| 47 | 47 |
CHARACTER(LCHARS) :: FileSpline,TmpChar |
| 48 | 48 |
|
| 49 | 49 |
|
| 50 |
! LOGICAL :: FAlign=.FALSE., FRot=.FALSE.
|
|
| 50 |
! LOGICAL :: FAlign=.FALSE., FRot=.FALSE.
|
|
| 51 | 51 |
|
| 52 | 52 |
! We will calculate the length of the path, in MW coordinates... |
| 53 | 53 |
! this is done is a stupid way: we interpolate the zmatrix values, |
| src/Hinvup_BFGS.f90 (revision 4) | ||
|---|---|---|
| 15 | 15 |
REAL(KREAL), INTENT(IN) :: ds(Nfree) |
| 16 | 16 |
REAL(KREAL), INTENT(IN) :: dGrad(Nfree) |
| 17 | 17 |
|
| 18 |
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|
| 19 |
! |
|
| 20 |
! This subroutine updates the inverse of the Hessian |
|
| 21 |
! using a BFGS formula. |
|
| 22 |
! |
|
| 23 |
! input: |
|
| 24 |
! - nfree: number of coordinates |
|
| 25 |
! - ds: GeomOld-Geom (expressed in the optimization coord.) |
|
| 26 |
! - dgrad: GradOld-Grad (expressed in the optimization coord.) |
|
| 27 |
! |
|
| 28 |
! input/output: |
|
| 29 |
! - Hess: old Hessian in input, replaced by updated Hessian on output |
|
| 30 |
! |
|
| 31 |
! Notations for the formula: |
|
| 32 |
! s=Geom-GeomOld |
|
| 33 |
! y=Grad-GradOld |
|
| 34 |
! W= (H)-1 ie the inverse matrix that we are updating |
|
| 35 |
! W+=W - [s(yT)W+Wy(sT)]/(y,s)+ {1+[(y,Wy)/(y,s)]}(s(sT))/(y,s)
|
|
| 36 |
! |
|
| 37 |
!!!!!!!!!!!!!!!! |
|
| 38 |
! |
|
| 39 |
! Due to the dual properties, this routine is also called to compute |
|
| 40 |
! the DFP update for the Hessian |
|
| 41 |
! |
|
| 42 |
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|
| 43 |
|
|
| 44 |
|
|
| 18 | 45 |
REAL(KREAL), ALLOCATABLE :: HessOld(:,:) , Wy(:) |
| 19 | 46 |
REAL(KREAL), ALLOCATABLE :: Hup(:,:),ssT(:,:) |
| 20 | 47 |
REAL(KREAL), ALLOCATABLE :: HupInt(:) |
| ... | ... | |
| 34 | 61 |
HupInt(nfree), & |
| 35 | 62 |
Num(nfree,nfree), Num2(nfree,nfree),Wy(Nfree),ssT(Nfree,nFree) ) |
| 36 | 63 |
|
| 37 |
! Notations |
|
| 38 |
! s=Geom-GeomOld |
|
| 39 |
! y=Grad-GradOld |
|
| 40 |
! W= (H)-1 ie the inverse matrix that we are updating |
|
| 41 |
! W+=W - [s(yT)W+Wy(sT)]/(y,s)+ {1+[(y,Wy)/(y,s)]}(s(sT))/(y,s)
|
|
| 42 | 64 |
|
| 43 | 65 |
|
| 44 | 66 |
HessOld=Hess |
| src/Step_DIIS_all.f90 (revision 4) | ||
|---|---|---|
| 3 | 3 |
!C Grad = input gradient vector. |
| 4 | 4 |
!C Geom_new = New Geometry. |
| 5 | 5 |
|
| 6 |
SUBROUTINE Step_diis_all(NGeomF,IGeom,Step,Geom,Grad,HP,HEAT,Hess,NCoord,allocation_flag,Tangent)
|
|
| 6 |
SUBROUTINE Step_diis_all(NGeomF,IGeom,Step,Geom,Grad,HP,HEAT,Hess,NCoord,allocation_flag,Tangent)
|
|
| 7 | 7 |
! IMPLICIT DOUBLE PRECISION (A-H,O-Z) |
| 8 | 8 |
|
| 9 |
USE Io_module, only : IoOut |
|
| 10 |
USE Path_module, only : Vfree |
|
| 11 |
|
|
| 12 |
IMPLICIT NONE |
|
| 13 |
INTEGER, parameter :: KINT = kind(1) |
|
| 14 |
INTEGER, parameter :: KREAL = kind(1.0d0) |
|
| 9 |
USE Io_module, only : IoOut |
|
| 10 |
USE Path_module, only : Vfree |
|
| 11 |
|
|
| 12 |
IMPLICIT NONE |
|
| 13 |
INTEGER, parameter :: KINT = kind(1) |
|
| 14 |
INTEGER, parameter :: KREAL = kind(1.0d0) |
|
| 15 |
|
|
| 16 |
! INCLUDE 'SIZES' |
|
| 17 |
|
|
| 18 |
INTEGER(KINT) :: NGeomF,IGeom |
|
| 19 |
INTEGER(KINT), INTENT(IN) :: NCoord |
|
| 20 |
|
|
| 21 |
REAL(KREAL) :: Geom_new(NCoord),Geom(NCoord),Grad(NCoord) |
|
| 22 |
REAL(KREAL) :: Hess(NCoord*NCoord),Step(NCoord) |
|
| 23 |
REAL(KREAL) :: HEAT,HP |
|
| 24 |
LOGICAL :: allocation_flag |
|
| 25 |
REAL(KREAL), INTENT(INOUT) :: Tangent(Ncoord) |
|
| 15 | 26 |
|
| 16 |
! INCLUDE 'SIZES' |
|
| 17 |
|
|
| 18 |
INTEGER(KINT) :: NGeomF,IGeom |
|
| 19 |
REAL(KREAL) :: Geom_new(NCoord),Geom(NCoord),Grad(NCoord) |
|
| 20 |
REAL(KREAL) :: Hess(NCoord*NCoord),Step(NCoord) |
|
| 21 |
REAL(KREAL) :: HEAT,HP |
|
| 22 |
LOGICAL :: allocation_flag |
|
| 23 |
INTEGER(KINT), INTENT(IN) :: NCoord |
|
| 24 |
REAL(KREAL), INTENT(INOUT) :: Tangent(Ncoord) |
|
| 25 |
|
|
| 26 | 27 |
!************************************************************************ |
| 27 | 28 |
!* * |
| 28 | 29 |
!* DIIS PERFORMS DIRECT INVERSION IN THE ITERATIVE SUBSPACE * |
| 29 | 30 |
!* * |
| 30 |
!* THIS INVOLVES SOLVING FOR C IN Geom(NEW) = Geom' - HG' *
|
|
| 31 |
!* THIS INVOLVES SOLVING FOR C IN Geom(NEW) = Geom' - HG' *
|
|
| 31 | 32 |
!* * |
| 32 | 33 |
!* WHERE Geom' = SUM(C(I)Geom(I), THE C COEFFICIENTES COMING FROM * |
| 33 | 34 |
!* * |
| ... | ... | |
| 35 | 36 |
!* | 1 0 | |-L | | 1 | * |
| 36 | 37 |
!* * |
| 37 | 38 |
!* WHERE B(I,J) =GRAD(I)H(T)HGRAD(J) GRAD(I) = GRADIENT ON CYCLE I * |
| 38 |
!* Hess = INVERSE HESSIAN *
|
|
| 39 |
!* Hess = INVERSE HESSIAN *
|
|
| 39 | 40 |
!* * |
| 40 | 41 |
!* REFERENCE * |
| 41 | 42 |
!* * |
| ... | ... | |
| 65 | 66 |
INTEGER(KINT), PARAMETER :: MRESET=15, M2=(MRESET+1)*(MRESET+1) !M2 = 256 |
| 66 | 67 |
REAL(KREAL), ALLOCATABLE, SAVE :: GeomSet(:,:),GradSet(:,:),ERR(:,:) ! MRESET*NCoord |
| 67 | 68 |
REAL(KREAL), SAVE :: ESET(MRESET) |
| 68 |
REAL(KREAL), ALLOCATABLE, SAVE :: DXTMP(:,:),GSAVE(:,:) !NCoord, why DXTMP has SAVE attribute??
|
|
| 69 |
REAL(KREAL), ALLOCATABLE, SAVE :: DXTMP(:,:),GSAVE(:,:) !NCoord, why DXTMP has SAVE attribute??
|
|
| 69 | 70 |
REAL(KREAL) :: B(M2),BS(M2),BST(M2) |
| 70 | 71 |
LOGICAL DEBUG, PRINT |
| 71 | 72 |
INTEGER(KINT), ALLOCATABLE, SAVE :: MSET(:) |
| 72 |
LOGICAL, ALLOCATABLE, SAVE :: FRST(:)
|
|
| 73 |
LOGICAL, ALLOCATABLE, SAVE :: FRST(:)
|
|
| 73 | 74 |
INTEGER(KINT) :: NDIIS, MPLUS, INV, ITERA, MM, NFree, I, J, K |
| 74 | 75 |
INTEGER(KINT) :: JJ, KJ, JNV, II, IONE, IJ, INK,ITmp, Isch, Idx |
| 75 | 76 |
REAL(KREAL) :: XMax, XNorm, S, DET, THRES, Norm |
| 76 |
REAL(KREAL), PARAMETER :: eps=1e-12
|
|
| 77 |
REAL(KREAL), PARAMETER :: crit=1e-8
|
|
| 78 |
REAL(KREAL), ALLOCATABLE :: Tanf(:) ! NCoord
|
|
| 79 |
REAL(KREAL), ALLOCATABLE :: HFree(:) ! NFree*NFree
|
|
| 80 |
REAL(KREAL), ALLOCATABLE :: Htmp(:,:) ! NCoord,NFree
|
|
| 81 |
REAL(KREAL), ALLOCATABLE :: Grad_free(:),Grad_new_free_inter(:),Step_free(:) ! NFree
|
|
| 82 |
REAL(KREAL), ALLOCATABLE :: Geom_free(:),Geom_new_free_inter(:),Geom_new_free(:) ! NFree
|
|
| 83 |
REAL(KREAL), ALLOCATABLE, SAVE :: GeomSet_free(:,:),GradSet_free(:,:)
|
|
| 84 |
|
|
| 77 |
REAL(KREAL), PARAMETER :: eps=1e-12
|
|
| 78 |
REAL(KREAL), PARAMETER :: crit=1e-8
|
|
| 79 |
REAL(KREAL), ALLOCATABLE :: Tanf(:) ! NCoord
|
|
| 80 |
REAL(KREAL), ALLOCATABLE :: HFree(:) ! NFree*NFree
|
|
| 81 |
REAL(KREAL), ALLOCATABLE :: Htmp(:,:) ! NCoord,NFree
|
|
| 82 |
REAL(KREAL), ALLOCATABLE :: Grad_free(:),Grad_new_free_inter(:),Step_free(:) ! NFree
|
|
| 83 |
REAL(KREAL), ALLOCATABLE :: Geom_free(:),Geom_new_free_inter(:),Geom_new_free(:) ! NFree
|
|
| 84 |
REAL(KREAL), ALLOCATABLE, SAVE :: GeomSet_free(:,:),GradSet_free(:,:)
|
|
| 85 |
|
|
| 85 | 86 |
DEBUG=.TRUE. |
| 86 | 87 |
PRINT=.TRUE. |
| 87 | 88 |
|
| 88 | 89 |
IF (PRINT) WRITE(*,'(/,'' BEGIN GDIIS '')') |
| 89 |
|
|
| 90 |
|
|
| 90 | 91 |
! Initialization |
| 91 | 92 |
IF (allocation_flag) THEN ! allocation_flag = .TRUE. at the begining and effective for all geometries in path. |
| 92 | 93 |
! FRST(IGeom) will be set to False in Space, so no need to modify it here |
| ... | ... | |
| 95 | 96 |
DEALLOCATE(GeomSet,GradSet,ERR,DXTMP,GSave,GeomSet_free,GradSet_free) |
| 96 | 97 |
RETURN |
| 97 | 98 |
ELSE |
| 98 |
! these allocated arrays need to be properly deallocated.
|
|
| 99 |
! these allocated arrays need to be properly deallocated.
|
|
| 99 | 100 |
IF (PRINT) WRITE(*,'(/,'' In FRST, GDIIS Alloc '')') |
| 100 | 101 |
ALLOCATE(GeomSet(NGeomF,MRESET*NCoord),GradSet(NGeomF,MRESET*NCoord),ERR(NGeomF,MRESET*NCoord)) |
| 101 |
ALLOCATE(GeomSet_free(NGeomF,MRESET*NCoord),GradSet_free(NGeomF,MRESET*NCoord))
|
|
| 102 |
ALLOCATE(GeomSet_free(NGeomF,MRESET*NCoord),GradSet_free(NGeomF,MRESET*NCoord))
|
|
| 102 | 103 |
ALLOCATE(DXTMP(NGeomF,NCoord),GSAVE(NGeomF,NCoord),MSET(NGeomF),FRST(NGeomF)) |
| 103 |
DO I=1,NGeomF
|
|
| 104 |
FRST(I) = .TRUE.
|
|
| 105 |
GeomSet(I,:) = 0.d0
|
|
| 106 |
GradSet(I,:) = 0.d0
|
|
| 107 |
ERR(I,:) = 0.d0
|
|
| 108 |
GeomSet_free(I,:) = 0.d0
|
|
| 109 |
GradSet_free(I,:) = 0.d0
|
|
| 110 |
DXTMP(I,:)=0.d0
|
|
| 111 |
GSAVE(I,:)=0.d0
|
|
| 112 |
END DO
|
|
| 113 |
MSET(:)=0
|
|
| 104 |
DO I=1,NGeomF
|
|
| 105 |
FRST(I) = .TRUE.
|
|
| 106 |
GeomSet(I,:) = 0.d0
|
|
| 107 |
GradSet(I,:) = 0.d0
|
|
| 108 |
ERR(I,:) = 0.d0
|
|
| 109 |
GeomSet_free(I,:) = 0.d0
|
|
| 110 |
GradSet_free(I,:) = 0.d0
|
|
| 111 |
DXTMP(I,:)=0.d0
|
|
| 112 |
GSAVE(I,:)=0.d0
|
|
| 113 |
END DO
|
|
| 114 |
MSET(:)=0
|
|
| 114 | 115 |
END IF |
| 115 |
allocation_flag = .FALSE.
|
|
| 116 |
allocation_flag = .FALSE.
|
|
| 116 | 117 |
END IF |
| 117 |
|
|
| 118 |
|
|
| 118 | 119 |
! Addded from here: |
| 119 |
Call FreeMv(NCoord,Vfree) ! VFree(Ncoord,Ncoord)
|
|
| 120 |
Call FreeMv(NCoord,Vfree) ! VFree(Ncoord,Ncoord)
|
|
| 120 | 121 |
! we orthogonalize Vfree to the tangent vector of this geom only if Tangent/=0.d0 |
| 121 | 122 |
Norm=sqrt(dot_product(Tangent,Tangent)) |
| 122 | 123 |
IF (Norm.GT.eps) THEN |
| ... | ... | |
| 126 | 127 |
Tangent=Tangent/Norm |
| 127 | 128 |
|
| 128 | 129 |
! We convert Tangent into Vfree only displacements. This is useless for now (2007.Apr.23) |
| 129 |
! as Vfree=Id matrix but it will be usefull as soon as we introduce constraints.
|
|
| 130 |
! as Vfree=Id matrix but it will be usefull as soon as we introduce constraints.
|
|
| 130 | 131 |
DO I=1,NCoord |
| 131 | 132 |
Tanf(I)=dot_product(reshape(Vfree(:,I),(/NCoord/)),Tangent) |
| 132 | 133 |
END DO |
| ... | ... | |
| 138 | 139 |
DO I=1,NCoord |
| 139 | 140 |
Norm=dot_product(reshape(vfree(:,I),(/NCoord/)),Tangent) |
| 140 | 141 |
Vfree(:,I)=Vfree(:,I)-Norm*Tangent |
| 141 |
END DO
|
|
| 142 |
END DO
|
|
| 142 | 143 |
|
| 143 | 144 |
Idx=0. |
| 144 | 145 |
! Schmidt orthogonalization of the Vfree vectors |
| ... | ... | |
| 156 | 157 |
Vfree(:,Idx)=Vfree(:,I)/sqrt(Norm) |
| 157 | 158 |
END IF |
| 158 | 159 |
END DO |
| 159 |
|
|
| 160 |
|
|
| 160 | 161 |
Print *, 'Idx=', Idx |
| 161 | 162 |
IF (Idx/= NCoord-1) THEN |
| 162 | 163 |
WRITE(*,*) "Pb in orthogonalizing Vfree to tangent for geom",IGeom |
| ... | ... | |
| 176 | 177 |
! We now calculate the new step |
| 177 | 178 |
! we project the hessian onto the free vectors |
| 178 | 179 |
ALLOCATE(HFree(NFree*NFree),Htmp(NCoord,NFree),Grad_free(NFree),Grad_new_free_inter(NFree)) |
| 179 |
ALLOCATE(Geom_free(NFree),Step_free(NFree),Geom_new_free_inter(NFree),Geom_new_free(NFree))
|
|
| 180 |
ALLOCATE(Geom_free(NFree),Step_free(NFree),Geom_new_free_inter(NFree),Geom_new_free(NFree))
|
|
| 180 | 181 |
DO J=1,NFree |
| 181 | 182 |
DO I=1,NCoord |
| 182 | 183 |
Htmp(I,J)=0.d0 |
| ... | ... | |
| 187 | 188 |
END DO |
| 188 | 189 |
DO J=1,NFree |
| 189 | 190 |
DO I=1,NFree |
| 190 |
HFree(I+((J-1)*NFree))=0.d0
|
|
| 191 |
HFree(I+((J-1)*NFree))=0.d0
|
|
| 191 | 192 |
DO K=1,NCoord |
| 192 |
HFree(I+((J-1)*NFree))=HFree(I+((J-1)*NFree))+Vfree(K,I)*Htmp(K,J)
|
|
| 193 |
HFree(I+((J-1)*NFree))=HFree(I+((J-1)*NFree))+Vfree(K,I)*Htmp(K,J)
|
|
| 193 | 194 |
END DO |
| 194 | 195 |
END DO |
| 195 | 196 |
END DO |
| 196 | 197 |
|
| 197 | 198 |
DO I=1,NFree |
| 198 | 199 |
Grad_free(I)=dot_product(reshape(Vfree(:,I),(/NCoord/)),Grad) |
| 199 |
Geom_free(I)=dot_product(reshape(Vfree(:,I),(/NCoord/)),Geom)
|
|
| 200 |
Geom_free(I)=dot_product(reshape(Vfree(:,I),(/NCoord/)),Geom)
|
|
| 200 | 201 |
END DO |
| 201 | 202 |
!Added Ends here.*********************************************** |
| 202 |
|
|
| 203 |
|
|
| 203 | 204 |
!C SPACE SIMPLY LOADS THE CURRENT VALUES OF Geom AND GRAD INTO |
| 204 | 205 |
!C THE ARRAYS GeomSet AND GradSet |
| 205 | 206 |
!C HEAT is never used, not even in Space_all(...) |
| 206 | 207 |
|
| 207 | 208 |
CALL Space_all(NGeomF,IGeom,MRESET,MSET,Geom,Grad,HEAT,NCoord,GeomSet,GradSet,ESET,FRST) |
| 208 |
|
|
| 209 |
|
|
| 209 | 210 |
IF (PRINT) WRITE(*,'(/,'' GDIIS after Space '')') |
| 210 |
|
|
| 211 |
DO J=1,MSet(IGeom)
|
|
| 211 |
|
|
| 212 |
DO J=1,MSet(IGeom)
|
|
| 212 | 213 |
DO K=1,NFree |
| 213 | 214 |
GradSet_free(IGeom,((J-1)*NFree)+K)=dot_product(reshape(Vfree(:,K),(/NCoord/)),& |
| 214 |
GradSet(IGeom,((J-1)*NCoord)+1:((J-1)*NCoord)+NCoord))
|
|
| 215 |
GeomSet_free(IGeom,((J-1)*NFree)+K)=dot_product(reshape(Vfree(:,K),(/NCoord/)),&
|
|
| 216 |
GeomSet(IGeom,((J-1)*NCoord)+1:((J-1)*NCoord)+NCoord))
|
|
| 215 |
GradSet(IGeom,((J-1)*NCoord)+1:((J-1)*NCoord)+NCoord))
|
|
| 216 |
GeomSet_free(IGeom,((J-1)*NFree)+K)=dot_product(reshape(Vfree(:,K),(/NCoord/)),&
|
|
| 217 |
GeomSet(IGeom,((J-1)*NCoord)+1:((J-1)*NCoord)+NCoord))
|
|
| 217 | 218 |
END DO |
| 218 |
END DO
|
|
| 219 |
END DO
|
|
| 219 | 220 |
!C |
| 220 | 221 |
!C INITIALIZE SOME VARIABLES AND CONSTANTS |
| 221 | 222 |
!C |
| 222 | 223 |
NDIIS = MSET(IGeom) |
| 223 | 224 |
MPLUS = MSET(IGeom) + 1 |
| 224 |
MM = MPLUS * MPLUS
|
|
| 225 |
MM = MPLUS * MPLUS
|
|
| 225 | 226 |
!C |
| 226 | 227 |
!C COMPUTE THE APPROXIMATE ERROR VECTORS |
| 227 | 228 |
!C |
| ... | ... | |
| 255 | 256 |
JJ = JJ + 1 |
| 256 | 257 |
B(JJ)=0.D0 |
| 257 | 258 |
DO 50 K=1,NFree |
| 258 |
!Print *, 'B(',JJ,')=', B(JJ) + ERR(IGeom,INV+K) * ERR(IGeom,JNV+K)
|
|
| 259 |
!Print *, 'B(',JJ,')=', B(JJ) + ERR(IGeom,INV+K) * ERR(IGeom,JNV+K)
|
|
| 259 | 260 |
50 B(JJ) = B(JJ) + ERR(IGeom,INV+K) * ERR(IGeom,JNV+K) |
| 260 | 261 |
|
| 261 | 262 |
! The following shifting is required to correct indices of B_ij elements in the GDIIS matrix. |
| 262 |
! The correction is needed because the last coloumn of the matrix contains all 1 and one zero.
|
|
| 263 |
! The correction is needed because the last coloumn of the matrix contains all 1 and one zero.
|
|
| 263 | 264 |
DO 60 I=MSET(IGeom)-1,1,-1 |
| 264 | 265 |
DO 60 J=MSET(IGeom),1,-1 |
| 265 | 266 |
60 B(I*MSET(IGeom)+J+I) = B(I*MSET(IGeom)+J) |
| ... | ... | |
| 298 | 299 |
JNV = (J-1) * MPLUS |
| 299 | 300 |
IF (J.EQ.ITERA) B(INV + J) = 0.D0 |
| 300 | 301 |
B(JNV + I) = B(INV + J) |
| 301 |
!Print *,'B(JNV + I)=',B(JNV + I)
|
|
| 302 |
!Print *,'B(JNV + I)=',B(JNV + I)
|
|
| 302 | 303 |
120 CONTINUE |
| 303 | 304 |
B(IONE) = 1.0D0 |
| 304 | 305 |
130 CONTINUE |
| ... | ... | |
| 318 | 319 |
!C |
| 319 | 320 |
DO 160 I=1,MPLUS |
| 320 | 321 |
II = MPLUS * (I-1) + I |
| 321 |
!Print *, 'B(',II,')=', B(II)
|
|
| 322 |
!Print *, 'GSave(',IGeom,',',I,')=', 1.D0 / DSQRT(1.D-20+DABS(B(II)))
|
|
| 322 |
!Print *, 'B(',II,')=', B(II)
|
|
| 323 |
!Print *, 'GSave(',IGeom,',',I,')=', 1.D0 / DSQRT(1.D-20+DABS(B(II)))
|
|
| 323 | 324 |
160 GSAVE(IGeom,I) = 1.D0 / DSQRT(1.D-20+DABS(B(II))) |
| 324 | 325 |
|
| 325 | 326 |
GSAVE(IGeom,MPLUS) = 1.D0 |
| 326 | 327 |
!Print *, 'GSave(',IGeom,',',MPlus,')=1.D0'
|
| 327 |
|
|
| 328 |
|
|
| 328 | 329 |
DO 170 I=1,MPLUS |
| 329 | 330 |
DO 170 J=1,MPLUS |
| 330 | 331 |
IJ = MPLUS * (I-1) + J |
| ... | ... | |
| 347 | 348 |
DO 190 I=1,MPLUS |
| 348 | 349 |
DO 190 J=1,MPLUS |
| 349 | 350 |
IJ = MPLUS * (I-1) + J |
| 350 |
!Print *, 'B(',IJ,')=', B(IJ)
|
|
| 351 |
!Print *, 'GSAVE(',IGeom,',',I,')=', GSAVE(IGeom,I)
|
|
| 352 |
!Print *, 'GSAVE(',IGeom,',',J,')=', GSAVE(IGeom,J)
|
|
| 353 |
!Print *, 'B(',IJ,')=', B(IJ) * GSAVE(I) * GSAVE(J)
|
|
| 351 |
!Print *, 'B(',IJ,')=', B(IJ)
|
|
| 352 |
!Print *, 'GSAVE(',IGeom,',',I,')=', GSAVE(IGeom,I)
|
|
| 353 |
!Print *, 'GSAVE(',IGeom,',',J,')=', GSAVE(IGeom,J)
|
|
| 354 |
!Print *, 'B(',IJ,')=', B(IJ) * GSAVE(I) * GSAVE(J)
|
|
| 354 | 355 |
190 B(IJ) = B(IJ) * GSAVE(IGeom,I) * GSAVE(IGeom,J) |
| 355 | 356 |
!C |
| 356 | 357 |
!C COMPUTE THE INTERMEDIATE INTERPOLATED PARAMETER AND GRADIENT VECTORS |
| ... | ... | |
| 362 | 363 |
DO 200 I=1,MSET(IGeom) |
| 363 | 364 |
INK = (I-1) * NFree + K |
| 364 | 365 |
Geom_new_free_inter(K) = Geom_new_free_inter(K) + B(MPLUS*MSET(IGeom)+I) * GeomSet_free(IGeom,INK) |
| 365 |
!Print *, 'Geom_new_free_inter(',K,')=', Geom_new_free_inter(K)
|
|
| 366 |
!Print *, 'B(MPLUS*MSET(',IGeom,')+',I,')=', B(MPLUS*MSET(IGeom)+I)
|
|
| 367 |
!Print *, 'GeomSet_free(',IGeom,',',INK,')=', GeomSet_free(IGeom,INK)
|
|
| 366 |
!Print *, 'Geom_new_free_inter(',K,')=', Geom_new_free_inter(K)
|
|
| 367 |
!Print *, 'B(MPLUS*MSET(',IGeom,')+',I,')=', B(MPLUS*MSET(IGeom)+I)
|
|
| 368 |
!Print *, 'GeomSet_free(',IGeom,',',INK,')=', GeomSet_free(IGeom,INK)
|
|
| 368 | 369 |
200 Grad_new_free_inter(K) = Grad_new_free_inter(K) + B(MPLUS*MSET(IGeom)+I) * GradSet_free(IGeom,INK) |
| 369 | 370 |
HP=0.D0 |
| 370 | 371 |
DO 210 I=1,MSET(IGeom) |
| ... | ... | |
| 395 | 396 |
Geom_new_free_inter(K) = Geom_free(K) |
| 396 | 397 |
240 Grad_new_free_inter(K) = Grad_free(K) |
| 397 | 398 |
ENDIF ! matches IF (XNORM.GT.2.D0 .OR. DABS(DET).LT. THRES) THEN, L378 |
| 398 |
|
|
| 399 |
! q_{m+1} = q'_{m+1} - H^{-1}g'_{m+1}
|
|
| 400 |
Geom_new_free=0.d0
|
|
| 401 |
DO I = 1, NFree
|
|
| 402 |
DO J = 1, NFree
|
|
| 403 |
! If Hinv=.False., then we need to invert Hess
|
|
| 404 |
!Geom_new_free(:) = Geom_new_free(:) + HFree(:,I)*Grad_new_free_inter(I)
|
|
| 405 |
Geom_new_free(J) = Geom_new_free(J) + HFree(I+((J-1)*NFree))*Grad_new_free_inter(I)
|
|
| 406 |
END DO
|
|
| 407 |
END DO
|
|
| 408 |
Geom_new_free(:) = Geom_new_free_inter(:) - Geom_new_free(:)
|
|
| 409 |
|
|
| 410 |
Step_free = Geom_new_free - Geom_free
|
|
| 399 |
|
|
| 400 |
! q_{m+1} = q'_{m+1} - H^{-1}g'_{m+1}
|
|
| 401 |
Geom_new_free=0.d0
|
|
| 402 |
DO I = 1, NFree
|
|
| 403 |
DO J = 1, NFree
|
|
| 404 |
! If Hinv=.False., then we need to invert Hess
|
|
| 405 |
!Geom_new_free(:) = Geom_new_free(:) + HFree(:,I)*Grad_new_free_inter(I)
|
|
| 406 |
Geom_new_free(J) = Geom_new_free(J) + HFree(I+((J-1)*NFree))*Grad_new_free_inter(I)
|
|
| 407 |
END DO
|
|
| 408 |
END DO
|
|
| 409 |
Geom_new_free(:) = Geom_new_free_inter(:) - Geom_new_free(:)
|
|
| 410 |
|
|
| 411 |
Step_free = Geom_new_free - Geom_free
|
|
| 411 | 412 |
|
| 412 |
Step = 0.d0
|
|
| 413 |
Step = 0.d0
|
|
| 413 | 414 |
DO I=1,NFree |
| 414 |
Step = Step + Step_free(I)*Vfree(:,I)
|
|
| 415 |
Step = Step + Step_free(I)*Vfree(:,I)
|
|
| 415 | 416 |
END DO |
| 416 | 417 |
|
| 417 | 418 |
DEALLOCATE(Hfree,Htmp,Grad_free,Grad_new_free_inter,Step_free,Geom_free) |
| 418 | 419 |
DEALLOCATE(Geom_new_free_inter,Geom_new_free) |
| 419 |
|
|
| 420 |
IF (PRINT) WRITE(*,'(/,'' END GDIIS '',/)')
|
|
| 421 |
|
|
| 420 |
|
|
| 421 |
IF (PRINT) WRITE(*,'(/,'' END GDIIS '',/)')
|
|
| 422 |
|
|
| 422 | 423 |
END SUBROUTINE Step_diis_all |
| src/Makefile (revision 4) | ||
|---|---|---|
| 5 | 5 |
# You might also have to edit the location of |
| 6 | 6 |
# some libraries (like mkl for ifort) |
| 7 | 7 |
|
| 8 |
Machine=gfortran
|
|
| 8 |
Machine=arqg
|
|
| 9 | 9 |
|
| 10 | 10 |
########################################### |
| 11 | 11 |
# # |
| ... | ... | |
| 32 | 32 |
|
| 33 | 33 |
ifeq ($(Machine),gfortran) |
| 34 | 34 |
# Flags for gfortran |
| 35 |
COMP=gfortran -g -Wall -fbounds-check
|
|
| 35 |
COMP=gfortran -fbounds-check -pedantic -std=gnu
|
|
| 36 | 36 |
#COMP=g95 -fbounds-check -g -ftrace=full -ftrace=frame -save-temps |
| 37 | 37 |
F90=${COMP}
|
| 38 | 38 |
F77=${COMP}
|
| ... | ... | |
| 136 | 136 |
# Flags for arq Ifort |
| 137 | 137 |
COMP=ifort |
| 138 | 138 |
F90=${COMP} -g -check bounds -check format -check uninit -traceback
|
| 139 |
#F90=${COMP} -g -check all -traceback
|
|
| 139 | 140 |
F77=${F90}
|
| 140 | 141 |
#F90=${COMP}
|
| 141 | 142 |
LINK=${COMP} -lguide -lpthread -L/usr/lib/ \
|
| ... | ... | |
| 147 | 148 |
|
| 148 | 149 |
ifeq ($(Machine),arqg) |
| 149 | 150 |
# Flags for arq GFortran |
| 150 |
COMP=gfortran |
|
| 151 |
COMP=gfortran -fbacktrace -fbounds-check -pedantic -std=gnu
|
|
| 151 | 152 |
F90=${COMP}
|
| 152 | 153 |
F77=${F90}
|
| 153 | 154 |
#F90=${COMP}
|
| ... | ... | |
| 253 | 254 |
Hupdate_MS.f90 \ |
| 254 | 255 |
Hupdate_Bofill.f90 \ |
| 255 | 256 |
Hinvup_BFGS.f90 \ |
| 257 |
Hinvup_DFP.f90 \ |
|
| 256 | 258 |
CalcTangent.f90 \ |
| 257 | 259 |
Step_RFO_all.f90 \ |
| 258 | 260 |
Step_GEDIIS.f90 \ |
| ... | ... | |
| 348 | 350 |
all: path utils |
| 349 | 351 |
|
| 350 | 352 |
install: |
| 351 |
mkdir ${PWD}/../exe
|
|
| 353 |
mkdir -p ${PWD}/../exe
|
|
| 352 | 354 |
ln -s ${PWD}/Path.exe ../exe/.
|
| 353 | 355 |
ln -s ${PWD}/../utils/xyz2scan ../exe/.
|
| 354 | 356 |
ln -s ${PWD}/../utils/xyz2path ../exe/.
|
| ... | ... | |
| 369 | 371 |
@echo "Utilities have been created." |
| 370 | 372 |
@echo "Make sure that they are in your PATH environment" |
| 371 | 373 |
|
| 372 |
xyz2scan: ../utils/Xyz2Scan.f |
|
| 373 |
${F90} -o ../utils/xyz2scan ../utils/Xyz2Scan.f
|
|
| 374 |
xyz2scan: ../utils/Xyz2Scan.f90
|
|
| 375 |
${F90} -o ../utils/xyz2scan ../utils/Xyz2Scan.f90
|
|
| 374 | 376 |
|
| 375 |
xyz2path: ../utils/Xyz2Path.f ../utils/Xyz2Path.param |
|
| 376 |
${F90} -o ../utils/xyz2path ../utils/Xyz2Path.f
|
|
| 377 |
xyz2path: ../utils/Xyz2Path.f90 ../utils/Xyz2Path.param
|
|
| 378 |
${F90} -o ../utils/xyz2path ../utils/Xyz2Path.f90
|
|
| 377 | 379 |
|
| 378 | 380 |
tgz: ${SRC0} ${SRC} Makefile ${EXAMPLES} ${MODSRC} ${SRCLAPACKD} ${SRCLAPACKU} ${SRCBLAS}
|
| 379 |
tar -cvf Path_${Version}.tar ${SRC0} ${SRC} Makefile ${EXAMPLES} ${MODSRC}
|
|
| 380 |
gzip Path_${Version}.tar
|
|
| 381 |
mv Path_${Version}.tar.gz Path_${Version}.tgz
|
|
| 382 |
@echo "Path_${Version}.tgz has been created."
|
|
| 381 |
tar -cvf OpenPath_${Version}.tar ${SRC0} ${SRC} Makefile ${EXAMPLES} ${MODSRC}
|
|
| 382 |
gzip OpenPath_${Version}.tar
|
|
| 383 |
mv OpenPath_${Version}.tar.gz OpenPath_${Version}.tgz
|
|
| 384 |
@echo "OpenPath_${Version}.tgz has been created."
|
|
| 383 | 385 |
|
| 384 | 386 |
lapack: ${OBJLAPACK}
|
| 385 | 387 |
ar rcs ./lapack/lapack.a ${OBJLAPACK}
|
| src/Step_DIIS.f90 (revision 4) | ||
|---|---|---|
| 2 | 2 |
!C XPARAM = input parameter vector (Geometry). |
| 3 | 3 |
!C Grad = input gradient vector. |
| 4 | 4 |
SUBROUTINE Step_DIIS(XP,XPARAM,GP,GRAD,HP,HEAT,Hess,NVAR,FRST) |
| 5 |
!SUBROUTINE DIIS(XPARAM,STEP,GRAD,HP,HEAT,Hess,NVAR,FRST)
|
|
| 5 |
!SUBROUTINE DIIS(XPARAM,STEP,GRAD,HP,HEAT,Hess,NVAR,FRST)
|
|
| 6 | 6 |
! IMPLICIT DOUBLE PRECISION (A-H,O-Z) |
| 7 | 7 |
IMPLICIT NONE |
| 8 | 8 |
integer, parameter :: KINT = kind(1) |
| ... | ... | |
| 12 | 12 |
|
| 13 | 13 |
INTEGER(KINT) :: NVAR |
| 14 | 14 |
REAL(KREAL) :: XP(NVAR), XPARAM(NVAR), GP(NVAR), GRAD(NVAR), Hess(NVAR*NVAR) |
| 15 |
!REAL(KREAL) :: Hess_inv(NVAR,NVAR),XPARAM_old(NVAR),STEP(NVAR)
|
|
| 15 |
!REAL(KREAL) :: Hess_inv(NVAR,NVAR),XPARAM_old(NVAR),STEP(NVAR)
|
|
| 16 | 16 |
REAL(KREAL) :: HEAT, HP |
| 17 | 17 |
LOGICAL :: FRST |
| 18 |
|
|
| 18 |
|
|
| 19 | 19 |
!************************************************************************ |
| 20 | 20 |
!* * |
| 21 | 21 |
!* DIIS PERFORMS DIRECT INVERSION IN THE ITERATIVE SUBSPACE * |
| ... | ... | |
| 70 | 70 |
PRINT=.TRUE. |
| 71 | 71 |
|
| 72 | 72 |
IF (PRINT) WRITE(*,'(/,'' BEGIN GDIIS '')') |
| 73 |
|
|
| 74 |
!XPARAM_old = XPARAM
|
|
| 75 |
|
|
| 73 |
|
|
| 74 |
!XPARAM_old = XPARAM
|
|
| 75 |
|
|
| 76 | 76 |
! Initialization |
| 77 | 77 |
IF (FRST) THEN |
| 78 | 78 |
! FRST will be set to False in Space, so no need to modify it here |
| ... | ... | |
| 122 | 122 |
DO I=1,MM |
| 123 | 123 |
B(I) = 1.D0 |
| 124 | 124 |
END DO |
| 125 |
|
|
| 125 |
|
|
| 126 | 126 |
! B_ij calculations from <e_i|e_j> |
| 127 | 127 |
JJ=0 |
| 128 | 128 |
INV=-NVAR |
| ... | ... | |
| 137 | 137 |
50 B(JJ) = B(JJ) + ERR(INV+K) * ERR(JNV+K) |
| 138 | 138 |
|
| 139 | 139 |
! The following shifting is required to correct indices of B_ij elements in the GDIIS matrix. |
| 140 |
! The correction is needed because the last coloumn of the matrix contains all 1 and one zero.
|
|
| 140 |
! The correction is needed because the last coloumn of the matrix contains all 1 and one zero.
|
|
| 141 | 141 |
DO 60 I=MSET-1,1,-1 |
| 142 | 142 |
DO 60 J=MSET,1,-1 |
| 143 | 143 |
60 B(I*MSET+J+I) = B(I*MSET+J) |
| ... | ... | |
| 194 | 194 |
II = MPLUS * (I-1) + I |
| 195 | 195 |
GSAVE(I) = 1.D0 / DSQRT(1.D-20+DABS(B(II))) |
| 196 | 196 |
END DO |
| 197 |
|
|
| 197 |
|
|
| 198 | 198 |
GSAVE(MPLUS) = 1.D0 |
| 199 | 199 |
DO I=1,MPLUS |
| 200 | 200 |
DO J=1,MPLUS |
| ... | ... | |
| 202 | 202 |
B(IJ) = B(IJ) * GSAVE(I) * GSAVE(J) |
| 203 | 203 |
END DO |
| 204 | 204 |
END DO |
| 205 |
|
|
| 205 |
|
|
| 206 | 206 |
IF (DEBUG) THEN |
| 207 | 207 |
! OUTPUT SCALED GDIIS MATRIX: |
| 208 | 208 |
WRITE(*,'(/5X,'' GDIIS MATRIX (SCALED)'')') |
| ... | ... | |
| 229 | 229 |
GP(K) = 0.D0 |
| 230 | 230 |
DO I=1,MSET |
| 231 | 231 |
INK = (I-1) * NVAR + K |
| 232 |
!Print *, 'B(',MPLUS*MSET+I,')=', B(MPLUS*MSET+I)
|
|
| 232 |
!Print *, 'B(',MPLUS*MSET+I,')=', B(MPLUS*MSET+I)
|
|
| 233 | 233 |
XP(K) = XP(K) + B(MPLUS*MSET+I) * XSET(INK) |
| 234 | 234 |
GP(K) = GP(K) + B(MPLUS*MSET+I) * GSET(INK) |
| 235 | 235 |
END DO |
| 236 |
END DO
|
|
| 237 |
|
|
| 236 |
END DO
|
|
| 237 |
|
|
| 238 | 238 |
HP=0.D0 |
| 239 | 239 |
DO I=1,MSET |
| 240 | 240 |
HP=HP+B(MPLUS*MSET+I)*ESET(I) |
| 241 | 241 |
END DO |
| 242 |
|
|
| 242 |
|
|
| 243 | 243 |
DO K=1,NVAR |
| 244 | 244 |
DX(K) = XPARAM(K) - XP(K) |
| 245 | 245 |
END DO |
| ... | ... | |
| 259 | 259 |
WRITE(*,*) "THE DIIS STEP WILL BE REPEATED WITH A SMALLER SPACE" |
| 260 | 260 |
END IF |
| 261 | 261 |
DO K=1,MM |
| 262 |
B(K) = BS(K)
|
|
| 262 |
B(K) = BS(K)
|
|
| 263 | 263 |
END DO |
| 264 | 264 |
NDIIS = NDIIS - 1 |
| 265 | 265 |
IF (NDIIS .GT. 0) GO TO 80 |
| ... | ... | |
| 269 | 269 |
GP(K) = GRAD(K) |
| 270 | 270 |
END DO |
| 271 | 271 |
ENDIF ! matches IF (XNORM.GT.2.D0 .OR. DABS(DET).LT. THRES) THEN |
| 272 |
|
|
| 272 |
|
|
| 273 | 273 |
! q_{m+1} = q'_{m+1} - H^{-1}g'_{m+1}
|
| 274 |
! Hess is a symmetric matrix.
|
|
| 275 |
!Hess_inv = 1.d0 ! to be deleted.
|
|
| 276 |
!Call GenInv(NVAR,Reshape(Hess,(/NVAR,NVAR/)),Hess_inv,NVAR) ! Implemented in Mat_util.f90
|
|
| 277 |
! H^{-1}g'_{m+1}
|
|
| 278 |
!Print *, 'Hess_inv='
|
|
| 279 |
!Print *, Hess_inv
|
|
| 280 |
!XPARAM=0.d0
|
|
| 281 |
!DO I=1, NVAR
|
|
| 282 |
!XPARAM(:) = XPARAM(:) + Hess_inv(:,I)*GP(I)
|
|
| 283 |
!END DO
|
|
| 284 |
!XPARAM(:) = XP(:) - XPARAM(:) ! now XPARAM is a new geometry.
|
|
| 285 |
|
|
| 286 |
!STEP is the difference between the new and old geometry and thus "step":
|
|
| 274 |
! Hess is a symmetric matrix.
|
|
| 275 |
!Hess_inv = 1.d0 ! to be deleted.
|
|
| 276 |
!Call GenInv(NVAR,Reshape(Hess,(/NVAR,NVAR/)),Hess_inv,NVAR) ! Implemented in Mat_util.f90
|
|
| 277 |
! H^{-1}g'_{m+1}
|
|
| 278 |
!Print *, 'Hess_inv='
|
|
| 279 |
!Print *, Hess_inv
|
|
| 280 |
!XPARAM=0.d0
|
|
| 281 |
!DO I=1, NVAR
|
|
| 282 |
!XPARAM(:) = XPARAM(:) + Hess_inv(:,I)*GP(I)
|
|
| 283 |
!END DO
|
|
| 284 |
!XPARAM(:) = XP(:) - XPARAM(:) ! now XPARAM is a new geometry.
|
|
| 285 |
|
|
| 286 |
!STEP is the difference between the new and old geometry and thus "step":
|
|
| 287 | 287 |
!STEP = XPARAM - XPARAM_old |
| 288 | 288 |
|
| 289 | 289 |
IF (PRINT) WRITE(*,'(/,'' END GDIIS '',/)') |
| src/Step_GEDIIS_All.f90 (revision 4) | ||
|---|---|---|
| 1 | 1 |
! Geom = input parameter vector (Geometry), Grad = input gradient vector, HEAT is Energy(Geom) |
| 2 | 2 |
SUBROUTINE Step_GEDIIS_All(NGeomF,IGeom,Step,Geom,Grad,HEAT,Hess,NCoord,allocation_flag,Tangent) |
| 3 | 3 |
!SUBROUTINE Step_GEDIIS(Geom_new,Geom,Grad,HEAT,Hess,NCoord,FRST) |
| 4 |
use Io_module
|
|
| 5 |
use Path_module, only : Nom, Atome, OrderInv, indzmat, Pi, Nat, Vfree
|
|
| 4 |
use Io_module
|
|
| 5 |
use Path_module, only : Nom, Atome, OrderInv, indzmat, Pi, Nat, Vfree
|
|
| 6 | 6 |
IMPLICIT NONE |
| 7 | 7 |
|
| 8 | 8 |
INTEGER(KINT) :: NGeomF,IGeom |
| 9 | 9 |
INTEGER(KINT), INTENT(IN) :: NCoord |
| 10 | 10 |
REAL(KREAL) :: Geom(NCoord), Grad(NCoord), Hess(NCoord*NCoord), Step(NCoord), Geom_new(NCoord) |
| 11 |
REAL(KREAL) :: HEAT ! HEAT= Energy
|
|
| 12 |
LOGICAL :: allocation_flag
|
|
| 13 |
REAL(KREAL), INTENT(INOUT) :: Tangent(Ncoord)
|
|
| 14 |
|
|
| 11 |
REAL(KREAL) :: HEAT ! HEAT= Energy
|
|
| 12 |
LOGICAL :: allocation_flag
|
|
| 13 |
REAL(KREAL), INTENT(INOUT) :: Tangent(Ncoord)
|
|
| 14 |
|
|
| 15 | 15 |
! MRESET = maximum number of iterations. |
| 16 | 16 |
INTEGER(KINT), PARAMETER :: MRESET=15, M2=(MRESET+1)*(MRESET+1) !M2 = 256 |
| 17 | 17 |
REAL(KREAL), ALLOCATABLE, SAVE :: GeomSet(:,:), GradSet(:,:) ! NGeomF,MRESET*NCoord |
| 18 | 18 |
REAL(KREAL), ALLOCATABLE, SAVE :: GSAVE(:,:) !NGeomF,NCoord |
| 19 |
REAL(KREAL), ALLOCATABLE, SAVE :: ESET(:,:)
|
|
| 20 |
REAL(KREAL) :: ESET_tmp(MRESET), B(M2), BS(M2), BST(M2), B_tmp(M2) ! M2=256
|
|
| 19 |
REAL(KREAL), ALLOCATABLE, SAVE :: ESET(:,:)
|
|
| 20 |
REAL(KREAL) :: ESET_tmp(MRESET), B(M2), BS(M2), BST(M2), B_tmp(M2) ! M2=256
|
|
| 21 | 21 |
LOGICAL :: DEBUG, PRINT, ci_lt_zero |
| 22 | 22 |
INTEGER(KINT), ALLOCATABLE, SAVE :: MSET(:) ! mth Iteration |
| 23 |
LOGICAL, ALLOCATABLE, SAVE :: FRST(:)
|
|
| 24 |
REAL(KREAL) :: ci(MRESET), ci_tmp(MRESET) ! MRESET = maximum number of iterations.
|
|
| 23 |
LOGICAL, ALLOCATABLE, SAVE :: FRST(:)
|
|
| 24 |
REAL(KREAL) :: ci(MRESET), ci_tmp(MRESET) ! MRESET = maximum number of iterations.
|
|
| 25 | 25 |
INTEGER(KINT) :: NGEDIIS, MPLUS, INV, ITERA, MM, cis_zero, IXX, JXX, MSET_minus_cis_zero |
| 26 | 26 |
INTEGER(KINT) :: I,J,K, JJ, KJ, JNV, II, IONE, IJ, INK,ITmp, IX, JX, KX, MSET_C_cis_zero |
| 27 |
INTEGER(KINT) :: current_size_B_mat, MyPointer, Iat, Isch, NFree, Idx
|
|
| 27 |
INTEGER(KINT) :: current_size_B_mat, MyPointer, Iat, Isch, NFree, Idx
|
|
| 28 | 28 |
REAL(KREAL) :: XMax, XNorm, S, DET, THRES, tmp, ER_star, ER_star_tmp, Norm |
| 29 |
REAL(KREAL), PARAMETER :: eps=1e-12
|
|
| 29 |
REAL(KREAL), PARAMETER :: eps=1e-12
|
|
| 30 | 30 |
REAL(KREAL), PARAMETER :: crit=1e-8 |
| 31 |
REAL(KREAL), ALLOCATABLE :: Tanf(:) ! NCoord
|
|
| 32 |
REAL(KREAL), ALLOCATABLE :: HFree(:) ! NFree*NFree
|
|
| 33 |
REAL(KREAL), ALLOCATABLE :: Htmp(:,:) ! NCoord,NFree
|
|
| 34 |
REAL(KREAL), ALLOCATABLE :: Grad_free(:), Step_free(:) ! NFree
|
|
| 35 |
REAL(KREAL), ALLOCATABLE :: Geom_free(:), Geom_new_free(:) ! NFree
|
|
| 36 |
REAL(KREAL), ALLOCATABLE, SAVE :: GeomSet_free(:,:), GradSet_free(:,:)
|
|
| 31 |
REAL(KREAL), ALLOCATABLE :: Tanf(:) ! NCoord
|
|
| 32 |
REAL(KREAL), ALLOCATABLE :: HFree(:) ! NFree*NFree
|
|
| 33 |
REAL(KREAL), ALLOCATABLE :: Htmp(:,:) ! NCoord,NFree
|
|
| 34 |
REAL(KREAL), ALLOCATABLE :: Grad_free(:), Step_free(:) ! NFree
|
|
| 35 |
REAL(KREAL), ALLOCATABLE :: Geom_free(:), Geom_new_free(:) ! NFree
|
|
| 36 |
REAL(KREAL), ALLOCATABLE, SAVE :: GeomSet_free(:,:), GradSet_free(:,:)
|
|
| 37 | 37 |
|
| 38 | 38 |
DEBUG=.TRUE. |
| 39 | 39 |
PRINT=.FALSE. |
| 40 |
|
|
| 40 |
|
|
| 41 | 41 |
IF (PRINT) WRITE(*,'(/,'' BEGIN Step_GEDIIS_ALL '')') |
| 42 |
|
|
| 42 |
|
|
| 43 | 43 |
! Initialization |
| 44 | 44 |
IF (allocation_flag) THEN |
| 45 | 45 |
! allocation_flag will be set to False in SPACE_GEDIIS, so no need to modify it here |
| ... | ... | |
| 50 | 50 |
ELSE |
| 51 | 51 |
IF (PRINT) WRITE(*,'(/,'' In allocation_flag, GEDIIS_ALL Alloc '')') |
| 52 | 52 |
ALLOCATE(GeomSet(NGeomF,MRESET*NCoord),GradSet(NGeomF,MRESET*NCoord),GSAVE(NGeomF,NCoord)) |
| 53 |
ALLOCATE(GeomSet_free(NGeomF,MRESET*NCoord),GradSet_free(NGeomF,MRESET*NCoord))
|
|
| 54 |
ALLOCATE(MSET(NGeomF),FRST(NGeomF),ESET(NGeomF,MRESET))
|
|
| 55 |
DO I=1,NGeomF
|
|
| 56 |
FRST(I) = .TRUE.
|
|
| 57 |
GeomSet(I,:) = 0.d0
|
|
| 58 |
GradSet(I,:) = 0.d0
|
|
| 59 |
GSAVE(I,:)=0.d0
|
|
| 60 |
GeomSet_free(I,:) = 0.d0
|
|
| 61 |
GradSet_free(I,:) = 0.d0
|
|
| 62 |
END DO
|
|
| 63 |
MSET(:)=0
|
|
| 53 |
ALLOCATE(GeomSet_free(NGeomF,MRESET*NCoord),GradSet_free(NGeomF,MRESET*NCoord))
|
|
| 54 |
ALLOCATE(MSET(NGeomF),FRST(NGeomF),ESET(NGeomF,MRESET))
|
|
| 55 |
DO I=1,NGeomF
|
|
| 56 |
FRST(I) = .TRUE.
|
|
| 57 |
GeomSet(I,:) = 0.d0
|
|
| 58 |
GradSet(I,:) = 0.d0
|
|
| 59 |
GSAVE(I,:)=0.d0
|
|
| 60 |
GeomSet_free(I,:) = 0.d0
|
|
| 61 |
GradSet_free(I,:) = 0.d0
|
|
| 62 |
END DO
|
|
| 63 |
MSET(:)=0
|
|
| 64 | 64 |
END IF |
| 65 |
allocation_flag = .FALSE.
|
|
| 65 |
allocation_flag = .FALSE.
|
|
| 66 | 66 |
END IF ! IF (allocation_flag) THEN |
| 67 |
|
|
| 68 |
! ADDED FROM HERE:
|
|
| 69 |
Call FreeMv(NCoord,Vfree) ! VFree(Ncoord,Ncoord), as of now, an Identity matrix.
|
|
| 67 |
|
|
| 68 |
! ADDED FROM HERE:
|
|
| 69 |
Call FreeMv(NCoord,Vfree) ! VFree(Ncoord,Ncoord), as of now, an Identity matrix.
|
|
| 70 | 70 |
! we orthogonalize Vfree to the tangent vector of this geom only if Tangent/=0.d0 |
| 71 | 71 |
Norm=sqrt(dot_product(Tangent,Tangent)) |
| 72 | 72 |
IF (Norm.GT.eps) THEN |
| ... | ... | |
| 76 | 76 |
Tangent=Tangent/Norm |
| 77 | 77 |
|
| 78 | 78 |
! We convert Tangent into Vfree only displacements. This is useless for now (2007.Apr.23) |
| 79 |
! as Vfree=Id matrix but it will be usefull as soon as we introduce constraints.
|
|
| 79 |
! as Vfree=Id matrix but it will be usefull as soon as we introduce constraints.
|
|
| 80 | 80 |
DO I=1,NCoord |
| 81 | 81 |
Tanf(I)=dot_product(reshape(Vfree(:,I),(/NCoord/)),Tangent) |
| 82 | 82 |
END DO |
| ... | ... | |
| 88 | 88 |
DO I=1,NCoord |
| 89 | 89 |
Norm=dot_product(reshape(vfree(:,I),(/NCoord/)),Tangent) |
| 90 | 90 |
Vfree(:,I)=Vfree(:,I)-Norm*Tangent |
| 91 |
END DO
|
|
| 91 |
END DO
|
|
| 92 | 92 |
|
| 93 | 93 |
Idx=0. |
| 94 | 94 |
! Schmidt orthogonalization of the Vfree vectors |
| ... | ... | |
| 106 | 106 |
Vfree(:,Idx)=Vfree(:,I)/sqrt(Norm) |
| 107 | 107 |
END IF |
| 108 | 108 |
END DO |
| 109 |
|
|
| 109 |
|
|
| 110 | 110 |
IF (Idx/= NCoord-1) THEN |
| 111 | 111 |
WRITE(*,*) "Pb in orthogonalizing Vfree to tangent for geom",IGeom |
| 112 | 112 |
WRITE(IOOut,*) "Pb in orthogonalizing Vfree to tangent for geom",IGeom |
| ... | ... | |
| 125 | 125 |
! We now calculate the new step |
| 126 | 126 |
! we project the hessian onto the free vectors |
| 127 | 127 |
ALLOCATE(HFree(NFree*NFree),Htmp(NCoord,NFree),Grad_free(NFree)) |
| 128 |
ALLOCATE(Geom_free(NFree),Step_free(NFree),Geom_new_free(NFree))
|
|
| 128 |
ALLOCATE(Geom_free(NFree),Step_free(NFree),Geom_new_free(NFree))
|
|
| 129 | 129 |
DO J=1,NFree |
| 130 | 130 |
DO I=1,NCoord |
| 131 | 131 |
Htmp(I,J)=0.d0 |
| ... | ... | |
| 136 | 136 |
END DO |
| 137 | 137 |
DO J=1,NFree |
| 138 | 138 |
DO I=1,NFree |
| 139 |
HFree(I+((J-1)*NFree))=0.d0
|
|
| 139 |
HFree(I+((J-1)*NFree))=0.d0
|
|
| 140 | 140 |
DO K=1,NCoord |
| 141 |
HFree(I+((J-1)*NFree))=HFree(I+((J-1)*NFree))+Vfree(K,I)*Htmp(K,J)
|
|
| 141 |
HFree(I+((J-1)*NFree))=HFree(I+((J-1)*NFree))+Vfree(K,I)*Htmp(K,J)
|
|
| 142 | 142 |
END DO |
| 143 | 143 |
END DO |
| 144 | 144 |
END DO |
| 145 | 145 |
|
| 146 | 146 |
DO I=1,NFree |
| 147 | 147 |
Grad_free(I)=dot_product(reshape(Vfree(:,I),(/NCoord/)),Grad) |
| 148 |
Geom_free(I)=dot_product(reshape(Vfree(:,I),(/NCoord/)),Geom)
|
|
| 148 |
Geom_free(I)=dot_product(reshape(Vfree(:,I),(/NCoord/)),Geom)
|
|
| 149 | 149 |
END DO |
| 150 | 150 |
!ADDED ENDS HERE.*********************************************** |
| 151 |
|
|
| 151 |
|
|
| 152 | 152 |
! SPACE_GEDIIS SIMPLY LOADS THE CURRENT VALUES OF Geom AND Grad INTO THE ARRAYS GeomSet |
| 153 | 153 |
! AND GradSet, MSET is set to zero and then 1 in SPACE_GEDIIS_All at first iteration. |
| 154 | 154 |
CALL SPACE_GEDIIS_All(NGeomF,IGeom,MRESET,MSET,Geom,Grad,HEAT,NCoord,GeomSet,GradSet,ESET,FRST) |
| 155 | 155 |
|
| 156 | 156 |
IF (PRINT) WRITE(*,'(/,'' GEDIIS after SPACE_GEDIIS_ALL '')') |
| 157 |
|
|
| 158 |
DO J=1,MSet(IGeom)
|
|
| 157 |
|
|
| 158 |
DO J=1,MSet(IGeom)
|
|
| 159 | 159 |
DO K=1,NFree |
| 160 | 160 |
GradSet_free(IGeom,((J-1)*NFree)+K)=dot_product(reshape(Vfree(:,K),(/NCoord/)),& |
| 161 |
GradSet(IGeom,((J-1)*NCoord)+1:((J-1)*NCoord)+NCoord))
|
|
| 162 |
GeomSet_free(IGeom,((J-1)*NFree)+K)=dot_product(reshape(Vfree(:,K),(/NCoord/)),&
|
|
| 163 |
GeomSet(IGeom,((J-1)*NCoord)+1:((J-1)*NCoord)+NCoord))
|
|
| 161 |
GradSet(IGeom,((J-1)*NCoord)+1:((J-1)*NCoord)+NCoord))
|
|
| 162 |
GeomSet_free(IGeom,((J-1)*NFree)+K)=dot_product(reshape(Vfree(:,K),(/NCoord/)),&
|
|
| 163 |
GeomSet(IGeom,((J-1)*NCoord)+1:((J-1)*NCoord)+NCoord))
|
|
| 164 | 164 |
END DO |
| 165 |
END DO
|
|
| 165 |
END DO
|
|
| 166 | 166 |
|
| 167 | 167 |
! INITIALIZE SOME VARIABLES AND CONSTANTS: |
| 168 | 168 |
NGEDIIS = MSET(IGeom) !MSET=mth iteration |
| ... | ... | |
| 180 | 180 |
JNV=JNV+NFree |
| 181 | 181 |
JJ = JJ + 1 |
| 182 | 182 |
B(JJ)=0.D0 |
| 183 |
DO K=1, NFree
|
|
| 184 |
B(JJ) = B(JJ) + (((GradSet_free(IGeom,INV+K)-GradSet_free(IGeom,JNV+K))* &
|
|
| 185 |
(GeomSet_free(IGeom,INV+K)-GeomSet_free(IGeom,JNV+K)))/2.D0)
|
|
| 186 |
END DO
|
|
| 183 |
DO K=1, NFree
|
|
| 184 |
B(JJ) = B(JJ) + (((GradSet_free(IGeom,INV+K)-GradSet_free(IGeom,JNV+K))* &
|
|
| 185 |
(GeomSet_free(IGeom,INV+K)-GeomSet_free(IGeom,JNV+K)))/2.D0)
|
|
| 186 |
END DO
|
|
| 187 | 187 |
END DO |
| 188 | 188 |
END DO |
| 189 | 189 |
|
| 190 | 190 |
! The following shifting is required to correct indices of B_ij elements in the GEDIIS matrix. |
| 191 |
! The correction is needed because the last coloumn of the matrix contains all 1 and one zero.
|
|
| 191 |
! The correction is needed because the last coloumn of the matrix contains all 1 and one zero.
|
|
| 192 | 192 |
DO I=MSET(IGeom)-1,1,-1 |
| 193 | 193 |
DO J=MSET(IGeom),1,-1 |
| 194 | 194 |
B(I*MSET(IGeom)+J+I) = B(I*MSET(IGeom)+J) |
| 195 | 195 |
END DO |
| 196 |
END DO
|
|
| 196 |
END DO
|
|
| 197 | 197 |
|
| 198 | 198 |
! For the last row and last column of GEDIIS matrix: |
| 199 | 199 |
DO I=1,MPLUS |
| ... | ... | |
| 201 | 201 |
B(MPLUS*MSET(IGeom)+I) = 1.D0 |
| 202 | 202 |
END DO |
| 203 | 203 |
B(MM) = 0.D0 |
| 204 |
|
|
| 205 |
DO I=1, MPLUS
|
|
| 206 |
!WRITE(*,'(10(1X,F20.4))') B((I-1)*MPLUS+1:I*(MPLUS))
|
|
| 207 |
END DO
|
|
| 208 |
|
|
| 204 |
|
|
| 205 |
DO I=1, MPLUS
|
|
| 206 |
!WRITE(*,'(10(1X,F20.4))') B((I-1)*MPLUS+1:I*(MPLUS))
|
|
| 207 |
END DO
|
|
| 208 |
|
|
| 209 | 209 |
! ELIMINATE ERROR VECTORS WITH THE LARGEST NORM: |
| 210 | 210 |
80 CONTINUE |
| 211 | 211 |
DO I=1,MM !MM = (MSET(IGeom)+1) * (MSET(IGeom)+1) |
| 212 | 212 |
BS(I) = B(I) !just a copy of the original B (GEDIIS) matrix |
| 213 | 213 |
END DO |
| 214 |
|
|
| 214 |
|
|
| 215 | 215 |
IF (NGEDIIS .NE. MSET(IGeom)) THEN |
| 216 | 216 |
DO II=1,MSET(IGeom)-NGEDIIS |
| 217 | 217 |
XMAX = -1.D10 |
| ... | ... | |
| 228 | 228 |
IONE = INV + I |
| 229 | 229 |
ENDIF |
| 230 | 230 |
END DO |
| 231 |
|
|
| 231 |
|
|
| 232 | 232 |
DO I=1,MPLUS |
| 233 | 233 |
INV = (I-1) * MPLUS |
| 234 | 234 |
DO J=1,MPLUS |
| ... | ... | |
| 236 | 236 |
IF (J.EQ.ITERA) B(INV + J) = 0.D0 |
| 237 | 237 |
B(JNV + I) = B(INV + J) |
| 238 | 238 |
END DO |
| 239 |
END DO
|
|
| 239 |
END DO
|
|
| 240 | 240 |
B(IONE) = 1.0D0 |
| 241 | 241 |
END DO |
| 242 |
END IF ! matches IF (NGEDIIS .NE. MSET(IGeom)) THEN
|
|
| 243 |
|
|
| 242 |
END IF ! matches IF (NGEDIIS .NE. MSET(IGeom)) THEN
|
|
| 243 |
|
|
| 244 | 244 |
! SCALE GEDIIS MATRIX BEFORE INVERSION: |
| 245 | 245 |
DO I=1,MPLUS |
| 246 | 246 |
II = MPLUS * (I-1) + I ! B(II)=diagonal elements of B matrix |
| 247 | 247 |
GSAVE(IGeom,I) = 1.D0 / DSQRT(1.D-20+DABS(B(II))) |
| 248 |
!Print *, 'GSAVE(',IGeom,',',I,')=', GSAVE(IGeom,I)
|
|
| 248 |
!Print *, 'GSAVE(',IGeom,',',I,')=', GSAVE(IGeom,I)
|
|
| 249 | 249 |
END DO |
| 250 | 250 |
GSAVE(IGeom,MPLUS) = 1.D0 |
| 251 | 251 |
DO I=1,MPLUS |
| ... | ... | |
| 256 | 256 |
END DO |
| 257 | 257 |
|
| 258 | 258 |
! INVERT THE GEDIIS MATRIX B: |
| 259 |
DO I=1, MPLUS
|
|
| 260 |
!WRITE(*,'(10(1X,F20.4))') B((I-1)*MPLUS+1:I*(MPLUS))
|
|
| 261 |
END DO
|
|
| 262 |
|
|
| 259 |
DO I=1, MPLUS
|
|
| 260 |
!WRITE(*,'(10(1X,F20.4))') B((I-1)*MPLUS+1:I*(MPLUS))
|
|
| 261 |
END DO
|
|
| 262 |
|
|
| 263 | 263 |
CALL MINV(B,MPLUS,DET) ! matrix inversion. |
| 264 |
|
|
| 265 |
DO I=1, MPLUS
|
|
| 266 |
!WRITE(*,'(10(1X,F20.16))') B((I-1)*MPLUS+1:I*(MPLUS))
|
|
| 267 |
END DO
|
|
| 268 |
|
|
| 264 |
|
|
| 265 |
DO I=1, MPLUS
|
|
| 266 |
!WRITE(*,'(10(1X,F20.16))') B((I-1)*MPLUS+1:I*(MPLUS))
|
|
| 267 |
END DO
|
|
| 268 |
|
|
| 269 | 269 |
DO I=1,MPLUS |
| 270 | 270 |
DO J=1,MPLUS |
| 271 | 271 |
IJ = MPLUS * (I-1) + J |
| 272 | 272 |
B(IJ) = B(IJ) * GSAVE(IGeom,I) * GSAVE(IGeom,J) |
| 273 | 273 |
END DO |
| 274 | 274 |
END DO |
| 275 |
|
|
| 276 |
! COMPUTE THE NEW INTERPOLATED PARAMETER VECTOR (Geometry):
|
|
| 277 |
ci=0.d0
|
|
| 278 |
ci_tmp=0.d0
|
|
| 279 |
|
|
| 280 |
ci_lt_zero= .FALSE.
|
|
| 281 |
DO I=1, MSET(IGeom)
|
|
| 282 |
DO J=1, MSET(IGeom) ! B matrix is read column-wise
|
|
| 283 |
ci(I)=ci(I)+B((J-1)*(MPLUS)+I)*ESET(IGeom,J) !ESET is energy set.
|
|
| 284 |
END DO
|
|
| 285 |
ci(I)=ci(I)+B((MPLUS-1)*(MPLUS)+I)
|
|
| 286 |
!Print *, 'NO ci < 0 yet, c(',I,')=', ci(I)
|
|
| 287 |
IF((ci(I) .LT. 0.0D0) .OR. (ci(I) .GT. 1.0D0)) THEN
|
|
| 288 |
ci_lt_zero=.TRUE.
|
|
| 289 |
EXIT
|
|
| 290 |
END IF
|
|
| 275 |
|
|
| 276 |
! COMPUTE THE NEW INTERPOLATED PARAMETER VECTOR (Geometry):
|
|
| 277 |
ci=0.d0
|
|
| 278 |
ci_tmp=0.d0
|
|
| 279 |
|
|
| 280 |
ci_lt_zero= .FALSE.
|
|
| 281 |
DO I=1, MSET(IGeom)
|
|
| 282 |
DO J=1, MSET(IGeom) ! B matrix is read column-wise
|
|
| 283 |
ci(I)=ci(I)+B((J-1)*(MPLUS)+I)*ESET(IGeom,J) !ESET is energy set.
|
|
| 284 |
END DO
|
|
| 285 |
ci(I)=ci(I)+B((MPLUS-1)*(MPLUS)+I)
|
|
| 286 |
!Print *, 'NO ci < 0 yet, c(',I,')=', ci(I)
|
|
| 287 |
IF((ci(I) .LT. 0.0D0) .OR. (ci(I) .GT. 1.0D0)) THEN
|
|
| 288 |
ci_lt_zero=.TRUE.
|
|
| 289 |
EXIT
|
|
| 290 |
END IF
|
|
| 291 | 291 |
END DO !matches DO I=1, MSET(IGeom) |
| 292 |
|
|
| 293 |
IF (ci_lt_zero) Then
|
|
| 294 |
cis_zero = 0
|
|
| 292 |
|
|
| 293 |
IF (ci_lt_zero) Then
|
|
| 294 |
cis_zero = 0
|
|
| 295 | 295 |
ER_star = 0.D0 |
| 296 |
ER_star_tmp = 1e32
|
|
| 297 |
|
|
| 298 |
! B_ij calculations from <B_ij=(g_i-g_j)(R_i-R_j)>, Full B matrix created first and then rows and columns are removed.
|
|
| 296 |
ER_star_tmp = 1e32
|
|
| 297 |
|
|
| 298 |
! B_ij calculations from <B_ij=(g_i-g_j)(R_i-R_j)>, Full B matrix created first and then rows and columns are removed.
|
|
| 299 | 299 |
JJ=0 |
| 300 | 300 |
INV=-NFree |
| 301 | 301 |
DO IX=1,MSET(IGeom) |
| ... | ... | |
| 305 | 305 |
JNV=JNV+NFree |
| 306 | 306 |
JJ = JJ + 1 |
| 307 | 307 |
BST(JJ)=0.D0 |
| 308 |
DO KX=1, NFree
|
|
| 309 |
BST(JJ) = BST(JJ) + (((GradSet_free(IGeom,INV+KX)-GradSet_free(IGeom,JNV+KX))* &
|
|
| 310 |
(GeomSet_free(IGeom,INV+KX)-GeomSet_free(IGeom,JNV+KX)))/2.D0)
|
|
| 311 |
END DO
|
|
| 308 |
DO KX=1, NFree
|
|
| 309 |
BST(JJ) = BST(JJ) + (((GradSet_free(IGeom,INV+KX)-GradSet_free(IGeom,JNV+KX))* &
|
|
| 310 |
(GeomSet_free(IGeom,INV+KX)-GeomSet_free(IGeom,JNV+KX)))/2.D0)
|
|
| 311 |
END DO
|
|
| 312 | 312 |
END DO |
| 313 |
END DO
|
|
| 314 |
|
|
| 315 |
DO I=1, (2**MSET(IGeom))-2 ! all (2**MSET(IGeom))-2 combinations of cis, except the one where all cis are .GT. 0 and .LT. 1
|
|
| 316 |
!Print *, 'Entering into DO I=1, (2**MSET(IGeom))-2 loop, MSET(IGeom)=', MSET(IGeom), ', I=', I
|
|
| 317 |
ci(:)=1.D0
|
|
| 318 |
! find out which cis are zero in I:
|
|
| 319 |
DO IX=1, MSET(IGeom)
|
|
| 320 |
JJ=IAND(I, 2**(IX-1))
|
|
| 321 |
IF(JJ .EQ. 0) Then
|
|
| 322 |
ci(IX)=0.D0
|
|
| 323 |
END IF
|
|
| 324 |
END DO
|
|
| 325 |
|
|
| 326 |
ci_lt_zero = .FALSE.
|
|
| 327 |
! B_ij calculations from <B_ij=(g_i-g_j)(R_i-R_j)>, Full B matrix created first and then rows and columns are removed.
|
|
| 328 |
DO IX=1, MSET(IGeom)*MSET(IGeom)
|
|
| 313 |
END DO
|
|
| 314 |
|
|
| 315 |
DO I=1, (2**MSET(IGeom))-2 ! all (2**MSET(IGeom))-2 combinations of cis, except the one where all cis are .GT. 0 and .LT. 1
|
|
| 316 |
!Print *, 'Entering into DO I=1, (2**MSET(IGeom))-2 loop, MSET(IGeom)=', MSET(IGeom), ', I=', I
|
|
| 317 |
ci(:)=1.D0
|
|
| 318 |
! find out which cis are zero in I:
|
|
| 319 |
DO IX=1, MSET(IGeom)
|
|
| 320 |
JJ=IAND(I, 2**(IX-1))
|
|
| 321 |
IF(JJ .EQ. 0) Then
|
|
| 322 |
ci(IX)=0.D0
|
|
| 323 |
END IF
|
|
| 324 |
END DO
|
|
| 325 |
|
|
| 326 |
ci_lt_zero = .FALSE.
|
|
| 327 |
! B_ij calculations from <B_ij=(g_i-g_j)(R_i-R_j)>, Full B matrix created first and then rows and columns are removed.
|
|
| 328 |
DO IX=1, MSET(IGeom)*MSET(IGeom)
|
|
| 329 | 329 |
B(IX) = BST(IX) !just a copy of the original B (GEDIIS) matrix |
| 330 | 330 |
END DO |
| 331 |
|
|
| 331 |
|
|
| 332 | 332 |
! Removal of KXth row and KXth column in order to accomodate cI to be zero: |
| 333 |
current_size_B_mat=MSET(IGeom) |
|
| 334 |
cis_zero = 0 |
|
| 335 |
! The bits of I (index of the upper loop 'DO I=1, (2**MSET(IGeom))-2'), gives which cis are zero. |
|
| 336 |
DO KX=1, MSET(IGeom) ! searching for each bit of I (index of the upper loop 'DO I=1, (2**MSET(IGeom))-2') |
|
| 337 |
IF (ci(KX) .EQ. 0.D0) Then !remove KXth row and KXth column |
|
| 338 |
cis_zero = cis_zero + 1 |
|
| 339 |
|
|
| 340 |
! First row removal: (B matrix is read column-wise) |
|
| 341 |
JJ=0 |
|
| 333 |
current_size_B_mat=MSET(IGeom) |
|
Formats disponibles : Unified diff