root / src / Calc_baker.f90 @ 2
Historique | Voir | Annoter | Télécharger (40,61 ko)
| 1 | 1 | equemene | SUBROUTINE Calc_baker(atome,r_cov,fact,frozen,IGeomRef) |
|---|---|---|---|
| 2 | 1 | equemene | ! |
| 3 | 1 | equemene | ! This subroutine analyses a geometry to construct the baker |
| 4 | 1 | equemene | ! delocalized internal coordinates |
| 5 | 1 | equemene | ! v1.0 |
| 6 | 1 | equemene | ! We use only one geometry |
| 7 | 1 | equemene | ! |
| 8 | 1 | equemene | Use Path_module, only : max_Z,NMaxL,Nom,MaxFroz,NFroz,Pi,Massat,a0,BMat_BakerT, & |
| 9 | 1 | equemene | Nat,NCoord,XyzGeomI,NGeomI,NGeomF,UMat,NPrim,BTransInv, & |
| 10 | 1 | equemene | IntCoordI,Coordinate,CurrentCoord,ScanCoord,BprimT,BBT, & |
| 11 | 1 | equemene | BBT_inv,Xprimitive,XPrimitiveF,Symmetry_elimination,UMat_local, & |
| 12 | 1 | equemene | BTransInv_local, UMatF,BTransInvF, & |
| 13 | 1 | equemene | indzmat, dzdc,renum,order,OrderInv |
| 14 | 1 | equemene | ! BMat_BakerT(3*Nat,NCoord), NCoord=3*Nat or NFree=3*Nat-6-Symmetry_elimination |
| 15 | 1 | equemene | ! depending upon the coordinate choice. IntCoordI(NGeomI,NCoord) where |
| 16 | 1 | equemene | ! UMat(NPrim,NCoord), NCoord number of vectors in UMat matrix, i.e. NCoord |
| 17 | 1 | equemene | ! Baker coordinates. NPrim is the number of primitive internal coordinates. |
| 18 | 1 | equemene | |
| 19 | 1 | equemene | Use Io_module |
| 20 | 1 | equemene | IMPLICIT NONE |
| 21 | 1 | equemene | |
| 22 | 1 | equemene | ! IGeom: Geometry index, IGeomRef: Reference Geometry. |
| 23 | 1 | equemene | INTEGER(KINT), INTENT(IN) :: atome(Nat),IGeomRef |
| 24 | 1 | equemene | REAL(KREAL), INTENT(IN) :: fact |
| 25 | 1 | equemene | REAL(KREAL), INTENT(IN) :: r_cov(0:Max_Z) |
| 26 | 1 | equemene | ! Frozen contains the indices of frozen atoms |
| 27 | 1 | equemene | INTEGER(KINT), INTENT(IN) :: Frozen(*) |
| 28 | 1 | equemene | |
| 29 | 1 | equemene | INTEGER(KINT), ALLOCATABLE :: LIAISONS(:,:) ! (Nat,0:NMaxL) |
| 30 | 1 | equemene | REAL(KREAL), ALLOCATABLE :: Geom(:,:) !(3,Nat) |
| 31 | 1 | equemene | ! NPrim is the number of primitive coordinates and NCoord is the number |
| 32 | 1 | equemene | ! of internal coordinates. BMat is actually (NPrim,3*Nat). |
| 33 | 1 | equemene | REAL(KREAL), ALLOCATABLE :: GMat(:,:) !(NPrim,NPrim) |
| 34 | 1 | equemene | ! EigVec(..) contains ALL eigevectors of BMat times BprimT, NOT only Baker Coordinate vectors. |
| 35 | 1 | equemene | REAL(KREAL), ALLOCATABLE :: EigVec(:,:), EigVal(:) ! EigVec(NPrim,NPrim) |
| 36 | 1 | equemene | REAL(KREAL), ALLOCATABLE :: x(:), y(:), z(:) |
| 37 | 1 | equemene | REAL(KREAL), ALLOCATABLE :: x_refGeom(:), y_refGeom(:), z_refGeom(:) |
| 38 | 1 | equemene | !REAL(KREAL), ALLOCATABLE :: V(:,:) ! (3*Nat,3*Nat) |
| 39 | 1 | equemene | !REAL(KREAL), ALLOCATABLE :: P(:) ! Used in Baker Coordinates: Matrix Inversion |
| 40 | 1 | equemene | REAL(KREAL), ALLOCATABLE :: UMat_tmp(:,:) |
| 41 | 1 | equemene | INTEGER(KINT) :: NbBonds, NbAngles, NbDihedrals, IGeom |
| 42 | 1 | equemene | |
| 43 | 1 | equemene | real(KREAL) :: vx,vy,vz,dist, Norm |
| 44 | 1 | equemene | real(KREAL) :: vx1,vy1,vz1,norm1 |
| 45 | 1 | equemene | real(KREAL) :: vx2,vy2,vz2,norm2 |
| 46 | 1 | equemene | real(KREAL) :: vx3,vy3,vz3,norm3 |
| 47 | 1 | equemene | real(KREAL) :: vx4,vy4,vz4,norm4 |
| 48 | 1 | equemene | real(KREAL) :: vx5,vy5,vz5,norm5 |
| 49 | 1 | equemene | real(KREAL) :: val,val_d, Q, T |
| 50 | 1 | equemene | |
| 51 | 1 | equemene | INTEGER(KINT) :: I,J, n1,n2,n3,n4,IAt,IL,JL,IFrag,ITmp, K, KMax |
| 52 | 1 | equemene | INTEGER(KINT) :: I0, IOld, IAtTmp, Izm, JAt, Kat, Lat, L, NOUT |
| 53 | 1 | equemene | |
| 54 | 1 | equemene | REAL(KREAL) :: sAngleIatIKat, sAngleIIatLat |
| 55 | 1 | equemene | REAL(KREAL) :: DiheTmp |
| 56 | 1 | equemene | |
| 57 | 1 | equemene | LOGICAL :: debug, bond, AddPrimitiveCoord, FAIL |
| 58 | 1 | equemene | |
| 59 | 1 | equemene | INTERFACE |
| 60 | 1 | equemene | function valid(string) result (isValid) |
| 61 | 1 | equemene | CHARACTER(*), intent(in) :: string |
| 62 | 1 | equemene | logical :: isValid |
| 63 | 1 | equemene | END function VALID |
| 64 | 1 | equemene | |
| 65 | 1 | equemene | |
| 66 | 1 | equemene | FUNCTION angle(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2) |
| 67 | 1 | equemene | use Path_module, only : Pi,KINT, KREAL |
| 68 | 1 | equemene | real(KREAL) :: v1x,v1y,v1z,norm1 |
| 69 | 1 | equemene | real(KREAL) :: v2x,v2y,v2z,norm2 |
| 70 | 1 | equemene | real(KREAL) :: angle |
| 71 | 1 | equemene | END FUNCTION ANGLE |
| 72 | 1 | equemene | |
| 73 | 1 | equemene | |
| 74 | 1 | equemene | |
| 75 | 1 | equemene | FUNCTION angle_d(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2,v3x,v3y,v3z,norm3) |
| 76 | 1 | equemene | use Path_module, only : Pi,KINT, KREAL |
| 77 | 1 | equemene | real(KREAL) :: v1x,v1y,v1z,norm1 |
| 78 | 1 | equemene | real(KREAL) :: v2x,v2y,v2z,norm2 |
| 79 | 1 | equemene | real(KREAL) :: v3x,v3y,v3z,norm3 |
| 80 | 1 | equemene | real(KREAL) :: angle_d,ca,sa |
| 81 | 1 | equemene | END FUNCTION ANGLE_D |
| 82 | 1 | equemene | |
| 83 | 1 | equemene | |
| 84 | 1 | equemene | |
| 85 | 1 | equemene | |
| 86 | 1 | equemene | SUBROUTINE Calc_Xprim(nat,x,y,z,Coordinate,NPrim,XPrimitive,XPrimRef) |
| 87 | 1 | equemene | ! |
| 88 | 1 | equemene | ! This subroutine uses the description of a list of Coordinates |
| 89 | 1 | equemene | ! to compute the values of the coordinates for a given geometry. |
| 90 | 1 | equemene | ! |
| 91 | 1 | equemene | !!!!!!!!!! |
| 92 | 1 | equemene | ! Input: |
| 93 | 1 | equemene | ! Na: INTEGER, Number of atoms |
| 94 | 1 | equemene | ! x,y,z(Na): REAL, cartesian coordinates of the considered geometry |
| 95 | 1 | equemene | ! Coordinate (Pointer(ListCoord)): description of the wanted coordiantes |
| 96 | 1 | equemene | ! NPrim, INTEGER: Number of coordinates to compute |
| 97 | 1 | equemene | ! |
| 98 | 1 | equemene | ! Optional: XPrimRef(NPrim) REAL: array that contains coordinates values for |
| 99 | 1 | equemene | ! a former geometry. Useful for Dihedral angles evolution... |
| 100 | 1 | equemene | |
| 101 | 1 | equemene | !!!!!!!!!!! |
| 102 | 1 | equemene | ! Output: |
| 103 | 1 | equemene | ! XPrimimite(NPrim) REAL: array that will contain the values of the coordinates |
| 104 | 1 | equemene | ! |
| 105 | 1 | equemene | !!!!!!!!! |
| 106 | 1 | equemene | |
| 107 | 1 | equemene | Use VarTypes |
| 108 | 1 | equemene | Use Io_module |
| 109 | 1 | equemene | Use Path_module, only : pi |
| 110 | 1 | equemene | |
| 111 | 1 | equemene | IMPLICIT NONE |
| 112 | 1 | equemene | |
| 113 | 1 | equemene | Type (ListCoord), POINTER :: Coordinate |
| 114 | 1 | equemene | INTEGER(KINT), INTENT(IN) :: Nat,NPrim |
| 115 | 1 | equemene | REAL(KREAL), INTENT(IN) :: x(Nat), y(Nat), z(Nat) |
| 116 | 1 | equemene | REAL(KREAL), INTENT(IN), OPTIONAL :: XPrimRef(NPrim) |
| 117 | 1 | equemene | REAL(KREAL), INTENT(OUT) :: XPrimitive(NPrim) |
| 118 | 1 | equemene | |
| 119 | 1 | equemene | END SUBROUTINE CALC_XPRIM |
| 120 | 1 | equemene | END INTERFACE |
| 121 | 1 | equemene | |
| 122 | 1 | equemene | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
| 123 | 1 | equemene | ! |
| 124 | 1 | equemene | !!!!!!!! PFL Debug |
| 125 | 1 | equemene | ! |
| 126 | 1 | equemene | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
| 127 | 1 | equemene | REAL(KREAL), ALLOCATABLE :: V(:,:) ! (NCoord,NCoord) |
| 128 | 1 | equemene | REAL(KREAL) :: tol |
| 129 | 1 | equemene | INTEGER(KINT) :: I1,I2,I3,I4 |
| 130 | 1 | equemene | LOGICAL :: FALOBBT,FALOBPrimt,FAloBBTInv |
| 131 | 1 | equemene | LOGICAL :: debugPFL |
| 132 | 1 | equemene | REAL(KREAL), ALLOCATABLE :: val_zmat(:,:) |
| 133 | 1 | equemene | INTEGER(KINT), ALLOCATABLE :: indzmat0(:,:) |
| 134 | 1 | equemene | |
| 135 | 1 | equemene | |
| 136 | 1 | equemene | debug=valid("Calc_baker")
|
| 137 | 1 | equemene | debugPFL=valid("BakerPFL")
|
| 138 | 1 | equemene | if (debug) WRITE(*,*) "============= Entering Cal_baker ==============" |
| 139 | 1 | equemene | |
| 140 | 1 | equemene | IF (Nat.le.2) THEN |
| 141 | 1 | equemene | WRITE(*,*) "I do not work for less than 2 atoms :-p" |
| 142 | 1 | equemene | RETURN |
| 143 | 1 | equemene | END IF |
| 144 | 1 | equemene | |
| 145 | 1 | equemene | IF (debug) THEN |
| 146 | 1 | equemene | WRITE(*,*) "DBG Calc_baker, geom" |
| 147 | 1 | equemene | DO I=1,Nat |
| 148 | 1 | equemene | WRITE(*,'(1X,I5,":",I3,3(1X,F15.3))') I,atome(I),XyzGeomI(IGeomRef,1:3,I) |
| 149 | 1 | equemene | END DO |
| 150 | 1 | equemene | END IF |
| 151 | 1 | equemene | |
| 152 | 1 | equemene | |
| 153 | 1 | equemene | |
| 154 | 1 | equemene | |
| 155 | 1 | equemene | ! we construct the list of bonds, valence angles and dihedral angles using this connectivity |
| 156 | 1 | equemene | ALLOCATE(Coordinate) |
| 157 | 1 | equemene | NULLIFY(Coordinate%next) |
| 158 | 1 | equemene | NbBonds=0 |
| 159 | 1 | equemene | NbAngles=0 |
| 160 | 1 | equemene | NbDihedrals=0 |
| 161 | 1 | equemene | |
| 162 | 1 | equemene | CurrentCoord => Coordinate |
| 163 | 1 | equemene | |
| 164 | 1 | equemene | ALLOCATE(Liaisons(Nat,0:NMaxL),Geom(3,Nat),x(Nat),y(Nat),z(Nat)) |
| 165 | 1 | equemene | ALLOCATE(x_refGeom(Nat),y_refGeom(Nat),z_refGeom(Nat)) |
| 166 | 1 | equemene | |
| 167 | 1 | equemene | x_refGeom(1:Nat) = XyzGeomI(IGeomRef,1,1:Nat) |
| 168 | 1 | equemene | y_refGeom(1:Nat) = XyzGeomI(IGeomRef,2,1:Nat) |
| 169 | 1 | equemene | z_refGeom(1:Nat) = XyzGeomI(IGeomRef,3,1:Nat) |
| 170 | 1 | equemene | |
| 171 | 1 | equemene | !!!!!!!!!!!!!!!!!!!! |
| 172 | 1 | equemene | ! |
| 173 | 1 | equemene | ! PFL Debug |
| 174 | 1 | equemene | ! |
| 175 | 1 | equemene | !!!!!!!!!!!!!!!!!!!! |
| 176 | 1 | equemene | if (debugPFL) THEN |
| 177 | 1 | equemene | WRITE(*,*) "DBG PFL Calc_BAKER - Calculating Zmat" |
| 178 | 1 | equemene | if (.not.ALLOCATED(indzmat)) THEN |
| 179 | 1 | equemene | ALLOCATE(indzmat(Nat,5)) |
| 180 | 1 | equemene | END IF |
| 181 | 1 | equemene | IF (.NOT.ALLOCATED(DzDc)) THEN |
| 182 | 1 | equemene | ALLOCATE(DzDc(3,nat,3,Nat)) |
| 183 | 1 | equemene | END IF |
| 184 | 1 | equemene | IF (.NOT.ALLOCATED(Val_zmat)) THEN |
| 185 | 1 | equemene | ALLOCATE(val_zmat(nat,3)) |
| 186 | 1 | equemene | END IF |
| 187 | 1 | equemene | IF (.NOT.ALLOCATED(indzmat0)) THEN |
| 188 | 1 | equemene | ALLOCATE(indzmat0(nat,5)) |
| 189 | 1 | equemene | END IF |
| 190 | 1 | equemene | ! We construct a Zmat |
| 191 | 1 | equemene | IF (Frozen(1).NE.0) THEN |
| 192 | 1 | equemene | Renum=.TRUE. |
| 193 | 1 | equemene | !!! remplcaer indzmat par indzmat0 puis renumerotation |
| 194 | 1 | equemene | call Calc_zmat_const_frag(Nat,atome,IndZmat0,val_zmat, & |
| 195 | 1 | equemene | x_refGeom,y_refGeom,z_refGeom, & |
| 196 | 1 | equemene | r_cov,fact,frozen) |
| 197 | 1 | equemene | ELSE |
| 198 | 1 | equemene | !no zmat... we create our own |
| 199 | 1 | equemene | call Calc_zmat_frag(Nat,atome,IndZmat0,val_zmat, & |
| 200 | 1 | equemene | x_refGeom,y_refGeom,z_refGeom, & |
| 201 | 1 | equemene | r_cov,fact) |
| 202 | 1 | equemene | END IF |
| 203 | 1 | equemene | |
| 204 | 1 | equemene | DO I=1,Nat |
| 205 | 1 | equemene | Order(IndZmat0(I,1))=I |
| 206 | 1 | equemene | OrderInv(I)=IndZmat0(I,1) |
| 207 | 1 | equemene | END DO |
| 208 | 1 | equemene | IndZmat(1,1)=Order(IndZmat0(1,1)) |
| 209 | 1 | equemene | IndZmat(2,1)=Order(IndZmat0(2,1)) |
| 210 | 1 | equemene | IndZmat(2,2)=Order(IndZmat0(2,2)) |
| 211 | 1 | equemene | IndZmat(3,1)=Order(IndZmat0(3,1)) |
| 212 | 1 | equemene | IndZmat(3,2)=Order(IndZmat0(3,2)) |
| 213 | 1 | equemene | IndZmat(3,3)=Order(IndZmat0(3,3)) |
| 214 | 1 | equemene | DO I=4,Nat |
| 215 | 1 | equemene | DO J=1,4 |
| 216 | 1 | equemene | IndZmat(I,J)=Order(IndZmat0(I,J)) |
| 217 | 1 | equemene | END DO |
| 218 | 1 | equemene | end do |
| 219 | 1 | equemene | |
| 220 | 1 | equemene | |
| 221 | 1 | equemene | WRITE(*,*) "DBG PFL CalcBaker, IndZmat0" |
| 222 | 1 | equemene | DO I=1,Nat |
| 223 | 1 | equemene | WRITE(*,*) I, indZmat0(I,:) |
| 224 | 1 | equemene | end do |
| 225 | 1 | equemene | |
| 226 | 1 | equemene | WRITE(*,*) "DBG PFL CalcBaker, IndZmat" |
| 227 | 1 | equemene | DO I=1,Nat |
| 228 | 1 | equemene | WRITE(*,*) I, indZmat(I,:) |
| 229 | 1 | equemene | end do |
| 230 | 1 | equemene | |
| 231 | 1 | equemene | |
| 232 | 1 | equemene | DO I=1,Nat |
| 233 | 1 | equemene | I1=indzmat0(I,1) |
| 234 | 1 | equemene | I2=indzmat0(I,2) |
| 235 | 1 | equemene | I3=indzmat0(I,3) |
| 236 | 1 | equemene | I4=indzmat0(I,4) |
| 237 | 1 | equemene | ! Adding BOND between at1 and at2 |
| 238 | 1 | equemene | WRITE(*,*) "DBG Indzmat",I,I1,I2,I3,I4 |
| 239 | 1 | equemene | if (I2.NE.0) THEN |
| 240 | 1 | equemene | NbBonds=NbBonds+1 |
| 241 | 1 | equemene | ! values of the coordinate (bond) is calculated for the reference geometry |
| 242 | 1 | equemene | ! whereas the coordinates (bonds) (CurrentCoord%At1, CurrentCoord%At2, etc.) |
| 243 | 1 | equemene | ! are determined using all geometries. vecteur is defined in bib_oper.f |
| 244 | 1 | equemene | Call vecteur(I1,I2,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2) |
| 245 | 1 | equemene | ! CurrentCoord%value=Norm2 |
| 246 | 1 | equemene | WRITE(*,'(1X,A,I5,A,2(1X,F12.6))') "Bond",NbBonds," Norm2,val_zmat",Norm2,val_zmat(I,1) |
| 247 | 1 | equemene | CurrentCoord%value=val_zmat(I,1) |
| 248 | 1 | equemene | CurrentCoord%SignDihedral = 0 |
| 249 | 1 | equemene | CurrentCoord%At1=I1 |
| 250 | 1 | equemene | CurrentCoord%At2=I2 |
| 251 | 1 | equemene | CurrentCoord%Type="BOND" |
| 252 | 1 | equemene | IF (.NOT.ASSOCIATED(CurrentCoord%next)) THEN |
| 253 | 1 | equemene | ALLOCATE(CurrentCoord%next) |
| 254 | 1 | equemene | NULLIFY(CurrentCoord%next%next) |
| 255 | 1 | equemene | END IF |
| 256 | 1 | equemene | CurrentCoord => CurrentCoord%next |
| 257 | 1 | equemene | END IF !IE.NE.0 |
| 258 | 1 | equemene | if (I3.NE.0) THEN |
| 259 | 1 | equemene | |
| 260 | 1 | equemene | NbAngles=NbAngles+1 |
| 261 | 1 | equemene | ! values of the coordinate (bond) is calculated for the reference geometry |
| 262 | 1 | equemene | ! whereas the coordinates (bonds) (CurrentCoord%At1, CurrentCoord%At2, etc.) are |
| 263 | 1 | equemene | ! determined using all geometries. |
| 264 | 1 | equemene | Call vecteur(i2,I3,x_refGeom,y_refGeom,z_refGeom,vx1,vy1,vz1,Norm1) |
| 265 | 1 | equemene | Call vecteur(i2,I1,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2) |
| 266 | 1 | equemene | ! CurrentCoord%value=angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2)*Pi/180. |
| 267 | 1 | equemene | WRITE(*,'(1X,A,I5,A,2(1X,F12.6))') "Angle",NbAngles," angle,val_zmat", & |
| 268 | 1 | equemene | angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2),val_zmat(I,2) |
| 269 | 1 | equemene | CurrentCoord%value=val_zmat(I,2)*Pi/180. |
| 270 | 1 | equemene | CurrentCoord%SignDihedral = 0 |
| 271 | 1 | equemene | CurrentCoord%At1=I1 |
| 272 | 1 | equemene | CurrentCoord%At2=I2 |
| 273 | 1 | equemene | CurrentCoord%At3=I3 |
| 274 | 1 | equemene | CurrentCoord%Type="ANGLE" |
| 275 | 1 | equemene | IF (.NOT.ASSOCIATED(CurrentCoord%next)) THEN |
| 276 | 1 | equemene | ALLOCATE(CurrentCoord%next) |
| 277 | 1 | equemene | NULLIFY(CurrentCoord%next%next) |
| 278 | 1 | equemene | END IF |
| 279 | 1 | equemene | CurrentCoord => CurrentCoord%next |
| 280 | 1 | equemene | END IF ! matches IF (I3.NE.0) |
| 281 | 1 | equemene | if (I4.NE.0) THEN |
| 282 | 1 | equemene | NbDihedrals=NbDihedrals+1 |
| 283 | 1 | equemene | Call vecteur(I2,I1,x_refGeom,y_refGeom,z_refGeom,vx1,vy1,vz1,Norm1) |
| 284 | 1 | equemene | Call vecteur(I2,I3,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2) |
| 285 | 1 | equemene | Call vecteur(I3,I4,x_refGeom,y_refGeom,z_refGeom,vx3,vy3,vz3,Norm3) |
| 286 | 1 | equemene | CALL produit_vect(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2, & |
| 287 | 1 | equemene | vx4,vy4,vz4,norm4) |
| 288 | 1 | equemene | CALL produit_vect(vx3,vy3,vz3,norm3,vx2,vy2,vz2,norm2, & |
| 289 | 1 | equemene | vx5,vy5,vz5,norm5) |
| 290 | 1 | equemene | ! CurrentCoord%value=angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, & |
| 291 | 1 | equemene | ! vx2,vy2,vz2,norm2)*Pi/180. |
| 292 | 1 | equemene | WRITE(*,'(1X,A,I5,A,2(1X,F12.6))') "Dihedral",NbDihedrals, & |
| 293 | 1 | equemene | " angle_d,val_zmat", & |
| 294 | 1 | equemene | angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, & |
| 295 | 1 | equemene | vx2,vy2,vz2,norm2) & |
| 296 | 1 | equemene | ,val_zmat(I,3) |
| 297 | 1 | equemene | |
| 298 | 1 | equemene | CurrentCoord%value = val_zmat(I,3)*Pi/180. |
| 299 | 1 | equemene | CurrentCoord%SignDihedral = int(sign(1.D0,val_zmat(I,3))) |
| 300 | 1 | equemene | CurrentCoord%At1=I1 |
| 301 | 1 | equemene | CurrentCoord%At2=I2 |
| 302 | 1 | equemene | CurrentCoord%At3=I3 |
| 303 | 1 | equemene | CurrentCoord%At4=I4 |
| 304 | 1 | equemene | CurrentCoord%Type="DIHEDRAL" |
| 305 | 1 | equemene | IF (.NOT.ASSOCIATED(CurrentCoord%next)) THEN |
| 306 | 1 | equemene | ALLOCATE(CurrentCoord%next) |
| 307 | 1 | equemene | NULLIFY(CurrentCoord%next%next) |
| 308 | 1 | equemene | END IF |
| 309 | 1 | equemene | CurrentCoord => CurrentCoord%next |
| 310 | 1 | equemene | END IF ! matches IF (I4.NE.0) |
| 311 | 1 | equemene | END DO ! matches DO I=1,Nat |
| 312 | 1 | equemene | deallocate(indzmat0) |
| 313 | 1 | equemene | |
| 314 | 1 | equemene | ELSE !! if (debugPFL) |
| 315 | 1 | equemene | |
| 316 | 1 | equemene | DO IGeom=1, NGeomI |
| 317 | 1 | equemene | !IGeom=1 ! need when using only first geometry. |
| 318 | 1 | equemene | x(1:Nat) = XyzGeomI(IGeom,1,1:Nat) |
| 319 | 1 | equemene | y(1:Nat) = XyzGeomI(IGeom,2,1:Nat) |
| 320 | 1 | equemene | z(1:Nat) = XyzGeomI(IGeom,3,1:Nat) |
| 321 | 1 | equemene | |
| 322 | 1 | equemene | ! First step: we calculate the connectivity |
| 323 | 1 | equemene | Liaisons=0 |
| 324 | 1 | equemene | |
| 325 | 1 | equemene | Call CalcCnct(Nat,atome,x,y,z,LIAISONS,r_cov,fact) |
| 326 | 1 | equemene | |
| 327 | 1 | equemene | IF (debug) THEN |
| 328 | 1 | equemene | WRITE(*,*) 'Connectivity done' |
| 329 | 1 | equemene | DO I=1,Nat |
| 330 | 1 | equemene | WRITE(*,'(1X,I5,":",I3,"=",15I4)') I,(Liaisons(I,J),J=0,NMaxL) |
| 331 | 1 | equemene | WRITE(*,'(1X,I5,":",I3,"=",15I4)') I,(Liaisons(I,:)) |
| 332 | 1 | equemene | END DO |
| 333 | 1 | equemene | DO I=1,Nat |
| 334 | 1 | equemene | WRITE(*,'(1X,I5,":",I3," r_cov=",F15.6)') I,atome(I),r_cov(atome(I)) |
| 335 | 1 | equemene | END DO |
| 336 | 1 | equemene | END IF |
| 337 | 1 | equemene | |
| 338 | 1 | equemene | DO I=1,Nat |
| 339 | 1 | equemene | IF (Liaisons(I,0).NE.0) THEN |
| 340 | 1 | equemene | DO J=1,Liaisons(I,0) |
| 341 | 1 | equemene | ! We add the bond |
| 342 | 1 | equemene | Iat=Liaisons(I,J) |
| 343 | 1 | equemene | IF (Iat.GT.I) THEN |
| 344 | 1 | equemene | ! Checking the uniqueness of bond. |
| 345 | 1 | equemene | ! Duplicate bonds should not be added. |
| 346 | 1 | equemene | AddPrimitiveCoord=.True. |
| 347 | 1 | equemene | ScanCoord=>Coordinate |
| 348 | 1 | equemene | DO WHILE (Associated(ScanCoord%next)) |
| 349 | 1 | equemene | IF (ScanCoord%Type .EQ. 'BOND') THEN |
| 350 | 1 | equemene | IF (ScanCoord%At1 .EQ. Iat) THEN |
| 351 | 1 | equemene | IF (ScanCoord%At2 .EQ. I) THEN |
| 352 | 1 | equemene | AddPrimitiveCoord=.False. |
| 353 | 1 | equemene | END IF |
| 354 | 1 | equemene | END IF |
| 355 | 1 | equemene | END IF |
| 356 | 1 | equemene | ScanCoord => ScanCoord%next |
| 357 | 1 | equemene | END DO |
| 358 | 1 | equemene | |
| 359 | 1 | equemene | IF (AddPrimitiveCoord) THEN |
| 360 | 1 | equemene | NbBonds=NbBonds+1 |
| 361 | 1 | equemene | ! values of the coordinate (bond) is calculated for the reference geometry |
| 362 | 1 | equemene | ! whereas the coordinates (bonds) (CurrentCoord%At1, CurrentCoord%At2, etc.) |
| 363 | 1 | equemene | ! are determined using all geometries. vecteur is defined in bib_oper.f |
| 364 | 1 | equemene | Call vecteur(I,Iat,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2) |
| 365 | 1 | equemene | CurrentCoord%value=Norm2 |
| 366 | 1 | equemene | CurrentCoord%SignDihedral = 0 |
| 367 | 1 | equemene | CurrentCoord%At1=Iat |
| 368 | 1 | equemene | CurrentCoord%At2=I |
| 369 | 1 | equemene | CurrentCoord%Type="BOND" |
| 370 | 1 | equemene | IF (.NOT.ASSOCIATED(CurrentCoord%next)) THEN |
| 371 | 1 | equemene | ALLOCATE(CurrentCoord%next) |
| 372 | 1 | equemene | NULLIFY(CurrentCoord%next%next) |
| 373 | 1 | equemene | END IF |
| 374 | 1 | equemene | CurrentCoord => CurrentCoord%next |
| 375 | 1 | equemene | END IF ! matches IF (AddPrimitiveCoord) THEN |
| 376 | 1 | equemene | END IF ! matches IF (Iat.GT.I) THEN |
| 377 | 1 | equemene | |
| 378 | 1 | equemene | ! Call Annul(Liaisons,IAt,I,Nat,NMaxL) |
| 379 | 1 | equemene | ! Liaisons(Iat,0)=Liaisons(Iat,0)-1 |
| 380 | 1 | equemene | |
| 381 | 1 | equemene | ! We add the valence angles |
| 382 | 1 | equemene | DO K=J+1,Liaisons(I,0) |
| 383 | 1 | equemene | Kat=Liaisons(I,K) |
| 384 | 1 | equemene | IF (Kat.GT.I) THEN |
| 385 | 1 | equemene | ! Checking the uniqueness of valence angle. |
| 386 | 1 | equemene | ! Duplicate bonds should not be added. |
| 387 | 1 | equemene | AddPrimitiveCoord=.True. |
| 388 | 1 | equemene | ScanCoord=>Coordinate |
| 389 | 1 | equemene | DO WHILE (Associated(ScanCoord%next)) |
| 390 | 1 | equemene | IF (ScanCoord%Type .EQ. 'ANGLE') THEN |
| 391 | 1 | equemene | IF (ScanCoord%At1 .EQ. Iat) THEN |
| 392 | 1 | equemene | IF (ScanCoord%At2 .EQ. I) THEN |
| 393 | 1 | equemene | IF (ScanCoord%At3 .EQ. Kat) THEN |
| 394 | 1 | equemene | AddPrimitiveCoord=.False. |
| 395 | 1 | equemene | END IF |
| 396 | 1 | equemene | END IF |
| 397 | 1 | equemene | END IF |
| 398 | 1 | equemene | END IF |
| 399 | 1 | equemene | ScanCoord => ScanCoord%next |
| 400 | 1 | equemene | END DO ! matches DO WHILE (Associated(ScanCoord%next)) |
| 401 | 1 | equemene | |
| 402 | 1 | equemene | IF (AddPrimitiveCoord) THEN |
| 403 | 1 | equemene | IF (debug) WRITE(*,*) "Adding angle:",Iat,I,Kat |
| 404 | 1 | equemene | NbAngles=NbAngles+1 |
| 405 | 1 | equemene | ! values of the coordinate (bond) is calculated for the reference geometry |
| 406 | 1 | equemene | ! whereas the coordinates (bonds) (CurrentCoord%At1, CurrentCoord%At2, etc.) are |
| 407 | 1 | equemene | ! determined using all geometries. |
| 408 | 1 | equemene | Call vecteur(i,kat,x_refGeom,y_refGeom,z_refGeom,vx1,vy1,vz1,Norm1) |
| 409 | 1 | equemene | Call vecteur(i,Iat,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2) |
| 410 | 1 | equemene | CurrentCoord%value=angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2)*Pi/180. |
| 411 | 1 | equemene | CurrentCoord%SignDihedral = 0 |
| 412 | 1 | equemene | sAngleIatIKat=abs(sin(angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2)*Pi/180.)) |
| 413 | 1 | equemene | CurrentCoord%At1=Iat |
| 414 | 1 | equemene | CurrentCoord%At2=I |
| 415 | 1 | equemene | CurrentCoord%At3=KAt |
| 416 | 1 | equemene | CurrentCoord%Type="ANGLE" |
| 417 | 1 | equemene | IF (.NOT.ASSOCIATED(CurrentCoord%next)) THEN |
| 418 | 1 | equemene | ALLOCATE(CurrentCoord%next) |
| 419 | 1 | equemene | NULLIFY(CurrentCoord%next%next) |
| 420 | 1 | equemene | END IF |
| 421 | 1 | equemene | CurrentCoord => CurrentCoord%next |
| 422 | 1 | equemene | END IF ! matches IF (AddPrimitiveCoord) THEN |
| 423 | 1 | equemene | ELSE ! matches IF (Kat.GT.I) THEN |
| 424 | 1 | equemene | ! Kat is less than I. If there is a bond between I and Kat, then this angle |
| 425 | 1 | equemene | ! has already been taken into account. |
| 426 | 1 | equemene | if (debug) WRITE(*,*) "Testing angle:",Iat,I,Kat |
| 427 | 1 | equemene | Bond=.FALSE. |
| 428 | 1 | equemene | DO L=1,Liaisons(Iat,0) |
| 429 | 1 | equemene | IF (Liaisons(Iat,L).EQ.Kat) Bond=.TRUE. |
| 430 | 1 | equemene | END DO |
| 431 | 1 | equemene | IF (.NOT.Bond) THEN |
| 432 | 1 | equemene | IF (debug) WRITE(*,*) "Not Bond Adding angle:",Iat,I,Kat |
| 433 | 1 | equemene | |
| 434 | 1 | equemene | ! Checking the uniqueness of valence angle. |
| 435 | 1 | equemene | ! Duplicate bonds should not be added. |
| 436 | 1 | equemene | AddPrimitiveCoord=.True. |
| 437 | 1 | equemene | ScanCoord=>Coordinate |
| 438 | 1 | equemene | DO WHILE (Associated(ScanCoord%next)) |
| 439 | 1 | equemene | IF (ScanCoord%Type .EQ. 'ANGLE') THEN |
| 440 | 1 | equemene | IF (ScanCoord%At1 .EQ. Iat) THEN |
| 441 | 1 | equemene | IF (ScanCoord%At2 .EQ. I) THEN |
| 442 | 1 | equemene | IF (ScanCoord%At3 .EQ. Kat) THEN |
| 443 | 1 | equemene | AddPrimitiveCoord=.False. |
| 444 | 1 | equemene | END IF |
| 445 | 1 | equemene | END IF |
| 446 | 1 | equemene | END IF |
| 447 | 1 | equemene | END IF |
| 448 | 1 | equemene | ScanCoord => ScanCoord%next |
| 449 | 1 | equemene | END DO |
| 450 | 1 | equemene | |
| 451 | 1 | equemene | IF (AddPrimitiveCoord) THEN |
| 452 | 1 | equemene | NbAngles=NbAngles+1 |
| 453 | 1 | equemene | ! values of the coordinate (bond) is calculated for the reference geometry |
| 454 | 1 | equemene | ! whereas the coordinates (bonds) (CurrentCoord%At1, CurrentCoord%At2, etc.) |
| 455 | 1 | equemene | ! are determined using all geometries. |
| 456 | 1 | equemene | Call vecteur(i,kat,x_refGeom,y_refGeom,z_refGeom,vx1,vy1,vz1,Norm1) |
| 457 | 1 | equemene | Call vecteur(i,Iat,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2) |
| 458 | 1 | equemene | CurrentCoord%value=angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2)*Pi/180. |
| 459 | 1 | equemene | CurrentCoord%SignDihedral = 0 |
| 460 | 1 | equemene | CurrentCoord%At1=Iat |
| 461 | 1 | equemene | CurrentCoord%At2=I |
| 462 | 1 | equemene | CurrentCoord%At3=KAt |
| 463 | 1 | equemene | CurrentCoord%Type="ANGLE" |
| 464 | 1 | equemene | IF (.NOT.ASSOCIATED(CurrentCoord%next)) THEN |
| 465 | 1 | equemene | ALLOCATE(CurrentCoord%next) |
| 466 | 1 | equemene | NULLIFY(CurrentCoord%next%next) |
| 467 | 1 | equemene | END IF |
| 468 | 1 | equemene | CurrentCoord => CurrentCoord%next |
| 469 | 1 | equemene | END IF ! matches IF (AddPrimitiveCoord) THEN |
| 470 | 1 | equemene | END IF ! matches IF (.NOT.Bond) THEN |
| 471 | 1 | equemene | END IF ! matches IF (Kat.GT.I) THEN, after else |
| 472 | 1 | equemene | END DO ! matches DO K=J+1,Liaisons(I,0) |
| 473 | 1 | equemene | END DO ! matches DO J=1,Liaisons(I,0) |
| 474 | 1 | equemene | END IF ! matches IF (Liaisons(I,0).NE.0) THEN |
| 475 | 1 | equemene | |
| 476 | 1 | equemene | ! We add dihedral angles. We do not concentrate on necessary angles, ie those that are |
| 477 | 1 | equemene | ! needed to build the molecule. Instead we collect all of them, and the choice of the best |
| 478 | 1 | equemene | ! ones will be done when diagonalizing the Baker matrix. |
| 479 | 1 | equemene | IF (Liaisons(I,0).GE.3) Then |
| 480 | 1 | equemene | DO J=1,Liaisons(I,0) |
| 481 | 1 | equemene | Iat=Liaisons(I,J) |
| 482 | 1 | equemene | ! values of the coordinate (bond) is calculated for the reference geometry |
| 483 | 1 | equemene | ! whereas the coordinates (bonds) (CurrentCoord%At1, CurrentCoord%At2, etc.) are |
| 484 | 1 | equemene | ! determined using all geometries. |
| 485 | 1 | equemene | Call vecteur(Iat,i,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2) |
| 486 | 1 | equemene | |
| 487 | 1 | equemene | DO K=J+1,Liaisons(I,0) |
| 488 | 1 | equemene | Kat=Liaisons(I,K) |
| 489 | 1 | equemene | Call vecteur(i,kat,x_refGeom,y_refGeom,z_refGeom,vx1,vy1,vz1,Norm1) |
| 490 | 1 | equemene | sAngleIatIKat=abs(sin(angle(vx1,vy1,vz1,Norm1,-vx2,-vy2,-vz2,Norm2)*Pi/180.)) |
| 491 | 1 | equemene | |
| 492 | 1 | equemene | DO L=K+1,Liaisons(I,0) |
| 493 | 1 | equemene | Lat=Liaisons(I,L) |
| 494 | 1 | equemene | if (debug) WRITE(*,*) "0-Dihe Kat,I,Iat:",Kat,I,Iat,angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2) |
| 495 | 1 | equemene | Call vecteur(iat,lat,x_refGeom,y_refGeom,z_refGeom,vx3,vy3,vz3,Norm3) |
| 496 | 1 | equemene | sAngleIIatLat=abs(sin(angle(vx3,vy3,vz3,Norm3,vx2,vy2,vz2,Norm2)*Pi/180.)) |
| 497 | 1 | equemene | if (debug) WRITE(*,*) "0-Dihe I,Iat,Lat:",angle(vx3,vy3,vz3,Norm3,vx2,vy2,vz2,Norm2) |
| 498 | 1 | equemene | IF ((sAngleIatIKat.ge.0.01d0).AND.(sAngleIIatLat.ge.0.01d0)) THEN |
| 499 | 1 | equemene | if (debug) WRITE(*,*) "Adding dihedral:",Kat,I,Iat,Lat |
| 500 | 1 | equemene | |
| 501 | 1 | equemene | ! Checking for the uniqueness of valence angle. |
| 502 | 1 | equemene | ! Duplicate bonds should not be added. |
| 503 | 1 | equemene | AddPrimitiveCoord=.True. |
| 504 | 1 | equemene | ScanCoord=>Coordinate |
| 505 | 1 | equemene | DO WHILE (Associated(ScanCoord%next)) |
| 506 | 1 | equemene | IF (ScanCoord%Type .EQ. 'DIHEDRAL') THEN |
| 507 | 1 | equemene | IF (ScanCoord%At1 .EQ. Kat) THEN |
| 508 | 1 | equemene | IF (ScanCoord%At2 .EQ. I) THEN |
| 509 | 1 | equemene | IF (ScanCoord%At3 .EQ. Iat) THEN |
| 510 | 1 | equemene | IF (ScanCoord%At4 .EQ. Lat) THEN |
| 511 | 1 | equemene | AddPrimitiveCoord=.False. |
| 512 | 1 | equemene | END IF |
| 513 | 1 | equemene | END IF |
| 514 | 1 | equemene | END IF |
| 515 | 1 | equemene | END IF |
| 516 | 1 | equemene | END IF |
| 517 | 1 | equemene | ScanCoord => ScanCoord%next |
| 518 | 1 | equemene | END DO ! matches DO WHILE (Associated(ScanCoord%next)) |
| 519 | 1 | equemene | |
| 520 | 1 | equemene | ! IF (AddPrimitiveCoord) THEN |
| 521 | 1 | equemene | ! NbDihedrals=NbDihedrals+1 |
| 522 | 1 | equemene | ! CALL produit_vect(vx3,vy3,vz3,norm3,vx2,vy2,vz2,norm2, & |
| 523 | 1 | equemene | ! vx5,vy5,vz5,norm5) |
| 524 | 1 | equemene | !CALL produit_vect(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2, & |
| 525 | 1 | equemene | ! vx4,vy4,vz4,norm4) |
| 526 | 1 | equemene | ! CurrentCoord%value=angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, & |
| 527 | 1 | equemene | ! vx2,vy2,vz2,norm2)*Pi/180. |
| 528 | 1 | equemene | ! CurrentCoord%SignDihedral = 0 |
| 529 | 1 | equemene | ! CurrentCoord%At1=Kat |
| 530 | 1 | equemene | ! CurrentCoord%At2=I |
| 531 | 1 | equemene | ! CurrentCoord%At3=IAt |
| 532 | 1 | equemene | ! CurrentCoord%At4=LAt |
| 533 | 1 | equemene | !WRITE(*,'(1X,":(dihed 1)",I5," - ",I5," - ",I5," - ",I5," = ",F15.6)') & |
| 534 | 1 | equemene | ! ScanCoord%At1,ScanCoord%At2,ScanCoord%At3,ScanCoord%At4,ScanCoord%Value*180./Pi |
| 535 | 1 | equemene | ! CurrentCoord%Type="DIHEDRAL" |
| 536 | 1 | equemene | ! IF (.NOT.ASSOCIATED(CurrentCoord%next)) THEN |
| 537 | 1 | equemene | ! ALLOCATE(CurrentCoord%next) |
| 538 | 1 | equemene | ! NULLIFY(CurrentCoord%next%next) |
| 539 | 1 | equemene | ! END IF |
| 540 | 1 | equemene | ! CurrentCoord => CurrentCoord%next |
| 541 | 1 | equemene | !END IF ! matches IF (AddPrimitiveCoord) THEN |
| 542 | 1 | equemene | |
| 543 | 1 | equemene | END IF ! matches IF ((sAngleIatIKat.ge.0.01d0).AND.(sAngleIIatLat.ge.0.01d0)) THEN |
| 544 | 1 | equemene | END DO ! matches DO L=K+1,Liaisons(I,0) |
| 545 | 1 | equemene | END DO ! matches DO K=J+1,Liaisons(I,0) |
| 546 | 1 | equemene | END DO ! matches DO J=1,Liaisons(I,0) |
| 547 | 1 | equemene | END IF !matches IF (Liaisons(I,0).GE.3) Then |
| 548 | 1 | equemene | |
| 549 | 1 | equemene | ! we now add other dihedrals |
| 550 | 1 | equemene | DO J=1,Liaisons(I,0) |
| 551 | 1 | equemene | Iat=Liaisons(I,J) |
| 552 | 1 | equemene | If (Iat.LE.I) Cycle |
| 553 | 1 | equemene | ! values of the coordinate (bond) is calculated for the reference geometry |
| 554 | 1 | equemene | ! whereas the coordinates (bonds) (CurrentCoord%At1, CurrentCoord%At2, etc.) are |
| 555 | 1 | equemene | ! determined using all geometries. |
| 556 | 1 | equemene | Call vecteur(Iat,i,x_refGeom,y_refGeom,z_refGeom,vx2,vy2,vz2,Norm2) |
| 557 | 1 | equemene | |
| 558 | 1 | equemene | DO K=1,Liaisons(I,0) |
| 559 | 1 | equemene | IF (Liaisons(I,K).EQ.Iat) Cycle |
| 560 | 1 | equemene | Kat=Liaisons(I,K) |
| 561 | 1 | equemene | Call vecteur(i,kat,x_refGeom,y_refGeom,z_refGeom,vx1,vy1,vz1,Norm1) |
| 562 | 1 | equemene | sAngleIatIKat=abs(sin(angle(vx1,vy1,vz1,Norm1,-vx2,-vy2,-vz2,Norm2)*Pi/180.)) |
| 563 | 1 | equemene | |
| 564 | 1 | equemene | DO L=1,Liaisons(Iat,0) |
| 565 | 1 | equemene | Lat=Liaisons(Iat,L) |
| 566 | 1 | equemene | IF (Lat.eq.I) Cycle |
| 567 | 1 | equemene | |
| 568 | 1 | equemene | if (debug) WRITE(*,*) "Dihe Kat,I,Iat:",angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2),sAngleIatIKat |
| 569 | 1 | equemene | Call vecteur(iat,lat,x_refGeom,y_refGeom,z_refGeom,vx,vy,vz,Norm) |
| 570 | 1 | equemene | |
| 571 | 1 | equemene | sAngleIIatLat=abs(sin(angle(vx,vy,vz,Norm,vx2,vy2,vz2,Norm2)*Pi/180.)) |
| 572 | 1 | equemene | if (debug) WRITE(*,*) "Dihe I,Iat,Lat:",angle(vx,vy,vz,Norm,vx2,vy2,vz2,Norm2) |
| 573 | 1 | equemene | ! we add the dihedral angle Kat-I-Iat-Lat |
| 574 | 1 | equemene | if (debug) WRITE(*,*) "Trying dihedral:",Kat,I,Iat,Lat,sAngleIatIKat,sAngleIIatLat |
| 575 | 1 | equemene | IF ((Lat.NE.Kat).AND.(Lat.Ne.I).AND.(sAngleIatIKat.ge.0.01d0) & |
| 576 | 1 | equemene | .AND.(sAngleIIatLat.ge.0.01d0)) THEN |
| 577 | 1 | equemene | if (debug) WRITE(*,*) "Adding dihedral:",Kat,I,Iat,Lat |
| 578 | 1 | equemene | |
| 579 | 1 | equemene | ! Checking for the uniqueness of valence angle. |
| 580 | 1 | equemene | ! Duplicate bonds should not be added. |
| 581 | 1 | equemene | AddPrimitiveCoord=.True. |
| 582 | 1 | equemene | ScanCoord=>Coordinate |
| 583 | 1 | equemene | DO WHILE (Associated(ScanCoord%next)) |
| 584 | 1 | equemene | IF (ScanCoord%Type .EQ. 'DIHEDRAL') THEN |
| 585 | 1 | equemene | IF (ScanCoord%At1 .EQ. Kat) THEN |
| 586 | 1 | equemene | IF (ScanCoord%At2 .EQ. I) THEN |
| 587 | 1 | equemene | IF (ScanCoord%At3 .EQ. Iat) THEN |
| 588 | 1 | equemene | IF (ScanCoord%At4 .EQ. Lat) THEN |
| 589 | 1 | equemene | AddPrimitiveCoord=.False. |
| 590 | 1 | equemene | END IF |
| 591 | 1 | equemene | END IF |
| 592 | 1 | equemene | END IF |
| 593 | 1 | equemene | END IF |
| 594 | 1 | equemene | END IF |
| 595 | 1 | equemene | ScanCoord => ScanCoord%next |
| 596 | 1 | equemene | END DO |
| 597 | 1 | equemene | |
| 598 | 1 | equemene | IF (AddPrimitiveCoord) THEN |
| 599 | 1 | equemene | NbDihedrals=NbDihedrals+1 |
| 600 | 1 | equemene | Call vecteur(iat,lat,x_refGeom,y_refGeom,z_refGeom,vx3,vy3,vz3,Norm3) |
| 601 | 1 | equemene | CALL produit_vect(vx3,vy3,vz3,norm3,vx2,vy2,vz2,norm2, & |
| 602 | 1 | equemene | vx5,vy5,vz5,norm5) |
| 603 | 1 | equemene | CALL produit_vect(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2, & |
| 604 | 1 | equemene | vx4,vy4,vz4,norm4) |
| 605 | 1 | equemene | CurrentCoord%value=angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, & |
| 606 | 1 | equemene | vx2,vy2,vz2,norm2)*Pi/180. |
| 607 | 1 | equemene | CurrentCoord%At1=Kat |
| 608 | 1 | equemene | CurrentCoord%At2=I |
| 609 | 1 | equemene | CurrentCoord%At3=IAt |
| 610 | 1 | equemene | CurrentCoord%At4=LAt |
| 611 | 1 | equemene | CurrentCoord%Type="DIHEDRAL" |
| 612 | 1 | equemene | CurrentCoord%SignDihedral=int(sign(1.d0,angle_d(vx4,vy4,vz4,norm4, & |
| 613 | 1 | equemene | vx5,vy5,vz5,norm5, vx2,vy2,vz2,norm2))) |
| 614 | 1 | equemene | WRITE(*,'(1X,":(dihed 2)",I5," - ",I5," - ",I5," - ",I5," = ",F15.6)') ScanCoord%At1,& |
| 615 | 1 | equemene | ScanCoord%At2, ScanCoord%At3,ScanCoord%At4,ScanCoord%Value*180./Pi |
| 616 | 1 | equemene | IF (.NOT.ASSOCIATED(CurrentCoord%next)) THEN |
| 617 | 1 | equemene | ALLOCATE(CurrentCoord%next) |
| 618 | 1 | equemene | NULLIFY(CurrentCoord%next%next) |
| 619 | 1 | equemene | END IF |
| 620 | 1 | equemene | CurrentCoord => CurrentCoord%next |
| 621 | 1 | equemene | END IF ! matches IF (AddPrimitiveCoord) THEN |
| 622 | 1 | equemene | END IF ! matches IF ((Lat.NE.Kat).AND.(Lat.Ne.I).AND.(sAngleIatIKat.ge.0.01d0) ..... |
| 623 | 1 | equemene | END DO ! matches DO L=1,Liaisons(Iat,0) |
| 624 | 1 | equemene | END DO ! matches DO K=1,Liaisons(I,0) |
| 625 | 1 | equemene | END DO ! matches DO J=1,Liaisons(I,0) |
| 626 | 1 | equemene | END DO ! end of loop over all atoms, DO I=1, Nat |
| 627 | 1 | equemene | END DO ! matches DO IGeom=1, NGeomI |
| 628 | 1 | equemene | |
| 629 | 1 | equemene | |
| 630 | 1 | equemene | |
| 631 | 1 | equemene | END IF ! if (debugPFL) |
| 632 | 1 | equemene | DEALLOCATE(Liaisons) |
| 633 | 1 | equemene | |
| 634 | 1 | equemene | WRITE(*,*) "I have found" |
| 635 | 1 | equemene | WRITE(*,*) NbBonds, " bonds" |
| 636 | 1 | equemene | WRITE(*,*) NbAngles, " angles" |
| 637 | 1 | equemene | WRITE(*,*) NbDihedrals, " dihedrals" |
| 638 | 1 | equemene | |
| 639 | 1 | equemene | WRITE(*,*) 'All in all...' |
| 640 | 1 | equemene | ScanCoord=>Coordinate |
| 641 | 1 | equemene | I=0 |
| 642 | 1 | equemene | DO WHILE (Associated(ScanCoord%next)) |
| 643 | 1 | equemene | I=I+1 |
| 644 | 1 | equemene | SELECT CASE (ScanCoord%Type) |
| 645 | 1 | equemene | CASE ('BOND')
|
| 646 | 1 | equemene | !WRITE(IOOUT,'(1X,I3,":",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1,ScanCoord%At2,ScanCoord%Value |
| 647 | 1 | equemene | WRITE(*,'(1X,I3,":(bond )",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1,ScanCoord%At2,ScanCoord%Value |
| 648 | 1 | equemene | CASE ('ANGLE')
|
| 649 | 1 | equemene | !WRITE(IOOUT,'(1X,I3,":",I5," - ",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1, & |
| 650 | 1 | equemene | ! ScanCoord%At2, ScanCoord%At3,ScanCoord%Value*180./Pi |
| 651 | 1 | equemene | WRITE(*,'(1X,I3,":(angle)",I5," - ",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1, & |
| 652 | 1 | equemene | ScanCoord%At2, ScanCoord%At3,ScanCoord%Value*180./Pi |
| 653 | 1 | equemene | CASE ('DIHEDRAL')
|
| 654 | 1 | equemene | !WRITE(IOOUT,'(1X,I3,":",I5," - ",I5," - ",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1,& |
| 655 | 1 | equemene | ! ScanCoord%At2, ScanCoord%At3,ScanCoord%At4,ScanCoord%Value*180./Pi |
| 656 | 1 | equemene | WRITE(*,'(1X,I3,":(dihed)",I5," - ",I5," - ",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1,& |
| 657 | 1 | equemene | ScanCoord%At2, ScanCoord%At3,ScanCoord%At4,ScanCoord%Value*180./Pi |
| 658 | 1 | equemene | END SELECT |
| 659 | 1 | equemene | ScanCoord => ScanCoord%next |
| 660 | 1 | equemene | END DO |
| 661 | 1 | equemene | |
| 662 | 1 | equemene | ! NPrim is the number of primitive coordinates. |
| 663 | 1 | equemene | NPrim=NbBonds+NbAngles+NbDihedrals |
| 664 | 1 | equemene | |
| 665 | 1 | equemene | ! Now calculating values of all primitive bonds for all geometries: |
| 666 | 1 | equemene | ALLOCATE(Xprimitive(NgeomI,NPrim),XPrimitiveF(NGeomF,NPrim)) |
| 667 | 1 | equemene | ALLOCATE(UMatF(NGeomF,NPrim,NCoord)) |
| 668 | 1 | equemene | ALLOCATE(BTransInvF(NGeomF,NCoord,3*Nat)) |
| 669 | 1 | equemene | |
| 670 | 1 | equemene | ! We calculate the Primitive values for the first geometry, this plays a special role |
| 671 | 1 | equemene | ! as it will serve as a reference for other geometries ! |
| 672 | 1 | equemene | |
| 673 | 1 | equemene | x(1:Nat) = XyzGeomI(1,1,1:Nat) |
| 674 | 1 | equemene | y(1:Nat) = XyzGeomI(1,2,1:Nat) |
| 675 | 1 | equemene | z(1:Nat) = XyzGeomI(1,3,1:Nat) |
| 676 | 1 | equemene | |
| 677 | 1 | equemene | Call Calc_Xprim(nat,x,y,z,Coordinate,NPrim,XPrimitive(1,:)) |
| 678 | 1 | equemene | |
| 679 | 1 | equemene | DO IGeom=2, NGeomI |
| 680 | 1 | equemene | |
| 681 | 1 | equemene | x(1:Nat) = XyzGeomI(IGeom,1,1:Nat) |
| 682 | 1 | equemene | y(1:Nat) = XyzGeomI(IGeom,2,1:Nat) |
| 683 | 1 | equemene | z(1:Nat) = XyzGeomI(IGeom,3,1:Nat) |
| 684 | 1 | equemene | |
| 685 | 1 | equemene | Call Calc_Xprim(nat,x,y,z,Coordinate,NPrim,XPrimitive(IGeom,:),Xprimitive(1,:)) |
| 686 | 1 | equemene | |
| 687 | 1 | equemene | ! ScanCoord=>Coordinate |
| 688 | 1 | equemene | ! I=0 ! index for the NPrim (NPrim is the number of primitive coordinates). |
| 689 | 1 | equemene | ! DO WHILE (Associated(ScanCoord%next)) |
| 690 | 1 | equemene | ! I=I+1 |
| 691 | 1 | equemene | ! SELECT CASE (ScanCoord%Type) |
| 692 | 1 | equemene | ! CASE ('BOND')
|
| 693 | 1 | equemene | ! Call vecteur(ScanCoord%At2,ScanCoord%At1,x,y,z,vx2,vy2,vz2,Norm2) |
| 694 | 1 | equemene | ! Xprimitive(IGeom,I) = Norm2 |
| 695 | 1 | equemene | ! CASE ('ANGLE')
|
| 696 | 1 | equemene | ! Call vecteur(ScanCoord%At2,ScanCoord%At3,x,y,z,vx1,vy1,vz1,Norm1) |
| 697 | 1 | equemene | ! Call vecteur(ScanCoord%At2,ScanCoord%At1,x,y,z,vx2,vy2,vz2,Norm2) |
| 698 | 1 | equemene | ! Xprimitive(IGeom,I) = angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2)*Pi/180. |
| 699 | 1 | equemene | ! CASE ('DIHEDRAL')
|
| 700 | 1 | equemene | |
| 701 | 1 | equemene | |
| 702 | 1 | equemene | ! Call vecteur(ScanCoord%At2,ScanCoord%At1,x,y,z,vx1,vy1,vz1,Norm1) |
| 703 | 1 | equemene | ! Call vecteur(ScanCoord%At2,ScanCoord%At3,x,y,z,vx2,vy2,vz2,Norm2) |
| 704 | 1 | equemene | ! Call vecteur(ScanCoord%At3,ScanCoord%At4,x,y,z,vx3,vy3,vz3,Norm3) |
| 705 | 1 | equemene | ! Call produit_vect(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2, & |
| 706 | 1 | equemene | ! vx4,vy4,vz4,norm4) |
| 707 | 1 | equemene | ! Call produit_vect(vx3,vy3,vz3,norm3,vx2,vy2,vz2,norm2, & |
| 708 | 1 | equemene | ! vx5,vy5,vz5,norm5) |
| 709 | 1 | equemene | |
| 710 | 1 | equemene | ! DiheTmp= angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, & |
| 711 | 1 | equemene | ! vx2,vy2,vz2,norm2) |
| 712 | 1 | equemene | ! Xprimitive(IGeom,I) = DiheTmp*Pi/180. |
| 713 | 1 | equemene | ! ! We treat large dihedral angles differently as +180=-180 mathematically and physically |
| 714 | 1 | equemene | ! ! but this causes lots of troubles when converting baker to cart. |
| 715 | 1 | equemene | ! ! So we ensure that large dihedral angles always have the same sign |
| 716 | 1 | equemene | ! if (abs(ScanCoord%SignDihedral).NE.1) THEN |
| 717 | 1 | equemene | ! ScanCoord%SignDihedral=Int(Sign(1.D0,DiheTmp)) |
| 718 | 1 | equemene | ! ELSE |
| 719 | 1 | equemene | ! If ((abs(DiheTmp).GE.170.D0).AND.(Sign(1.,DiheTmp)*ScanCoord%SignDihedral<0)) THEN |
| 720 | 1 | equemene | ! Xprimitive(IGeom,I) = DiheTmp*Pi/180.+ ScanCoord%SignDihedral*2.*Pi |
| 721 | 1 | equemene | ! END IF |
| 722 | 1 | equemene | ! END IF |
| 723 | 1 | equemene | ! ! Another test... will all this be consistent ??? |
| 724 | 1 | equemene | ! ! We want the shortest path, so we check that the change in dihedral angles is less than Pi: |
| 725 | 1 | equemene | ! IF (IGeom.GT.1) THEN |
| 726 | 1 | equemene | ! IF (abs(Xprimitive(IGeom,I)-XPrimitive(IGeom-1,I)).GE.Pi) THEN |
| 727 | 1 | equemene | ! if ((Xprimitive(IGeom,I)-XPrimitive(IGeom-1,I)).GE.Pi) THEN |
| 728 | 1 | equemene | ! Xprimitive(IGeom,I)=Xprimitive(IGeom,I)-2.d0*Pi |
| 729 | 1 | equemene | ! ELSE |
| 730 | 1 | equemene | ! Xprimitive(IGeom,I)=Xprimitive(IGeom,I)+2.d0*Pi |
| 731 | 1 | equemene | ! END IF |
| 732 | 1 | equemene | ! END IF |
| 733 | 1 | equemene | ! END IF |
| 734 | 1 | equemene | ! END SELECT |
| 735 | 1 | equemene | ! ScanCoord => ScanCoord%next |
| 736 | 1 | equemene | ! END DO ! matches DO WHILE (Associated(ScanCoord%next)) |
| 737 | 1 | equemene | END DO ! matches DO IGeom=1, NGeomI |
| 738 | 1 | equemene | |
| 739 | 1 | equemene | IF (debug) THEN |
| 740 | 1 | equemene | WRITE(*,*) "DBG Calc_Baker Xprimitives for all geometries" |
| 741 | 1 | equemene | DO IGeom=1, NGeomI |
| 742 | 1 | equemene | WRITE(*,*) "Geometry ",IGeom |
| 743 | 1 | equemene | ScanCoord=>Coordinate |
| 744 | 1 | equemene | I=0 ! index for the NPrim (NPrim is the number of primitive coordinates). |
| 745 | 1 | equemene | DO WHILE (Associated(ScanCoord%next)) |
| 746 | 1 | equemene | I=I+1 |
| 747 | 1 | equemene | SELECT CASE (ScanCoord%Type) |
| 748 | 1 | equemene | CASE ('BOND')
|
| 749 | 1 | equemene | WRITE(*,'(1X,I3,":",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1,ScanCoord%At2,Xprimitive(IGeom,I) |
| 750 | 1 | equemene | CASE ('ANGLE')
|
| 751 | 1 | equemene | WRITE(*,'(1X,I3,":",I5," - ",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1, & |
| 752 | 1 | equemene | ScanCoord%At2, ScanCoord%At3,Xprimitive(IGeom,I)*180./Pi |
| 753 | 1 | equemene | CASE ('DIHEDRAL')
|
| 754 | 1 | equemene | WRITE(*,'(1X,I3,":",I5," - ",I5," - ",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1,& |
| 755 | 1 | equemene | ScanCoord%At2, ScanCoord%At3,ScanCoord%At4,Xprimitive(IGeom,I)*180./Pi |
| 756 | 1 | equemene | END SELECT |
| 757 | 1 | equemene | ScanCoord => ScanCoord%next |
| 758 | 1 | equemene | END DO ! matches DO WHILE (Associated(ScanCoord%next)) |
| 759 | 1 | equemene | END DO ! matches DO IGeom=1, NGeomI |
| 760 | 1 | equemene | END IF |
| 761 | 1 | equemene | ! Here reference geometry is assigned to Geom. Earlier x,y,z(Nat) were assigned from XyzGeomI(...) |
| 762 | 1 | equemene | ! before the call of this (Calc_baker) subroutine and passed to this subroutine as arguments of the |
| 763 | 1 | equemene | ! subroutine. |
| 764 | 1 | equemene | ! PFL 19 Feb 2008 -> |
| 765 | 1 | equemene | ! In order to reproduce the B matrix in the Baker article, one has |
| 766 | 1 | equemene | ! to divide Geom by a0. |
| 767 | 1 | equemene | ! However, this is inconsitent with our way of storing geometries |
| 768 | 1 | equemene | ! so we do not do it ! |
| 769 | 1 | equemene | |
| 770 | 1 | equemene | Geom(1,:)=XyzGeomI(IGeomRef,1,1:Nat) ! XyzGeomI(NGeomI,3,Nat) |
| 771 | 1 | equemene | Geom(2,:)=XyzGeomI(IGeomRef,2,1:Nat) |
| 772 | 1 | equemene | Geom(3,:)=XyzGeomI(IGeomRef,3,1:Nat) |
| 773 | 1 | equemene | |
| 774 | 1 | equemene | ! <- PFL 19 Feb 2008 |
| 775 | 1 | equemene | |
| 776 | 1 | equemene | ! We now have to compute dq/dx... |
| 777 | 1 | equemene | ALLOCATE(BprimT(3*Nat,NPrim)) |
| 778 | 1 | equemene | BprimT=0.d0 |
| 779 | 1 | equemene | |
| 780 | 1 | equemene | ScanCoord=>Coordinate |
| 781 | 1 | equemene | I=0 |
| 782 | 1 | equemene | DO WHILE (Associated(ScanCoord%next)) |
| 783 | 1 | equemene | I=I+1 |
| 784 | 1 | equemene | SELECT CASE (ScanCoord%Type) |
| 785 | 1 | equemene | CASE ('BOND')
|
| 786 | 1 | equemene | CALL CONSTRAINTS_BONDLENGTH_DER(Nat,ScanCoord%at1,ScanCoord%AT2, & |
| 787 | 1 | equemene | Geom,BprimT(1,I)) |
| 788 | 1 | equemene | CASE ('ANGLE')
|
| 789 | 1 | equemene | CALL CONSTRAINTS_BONDANGLE_DER(Nat,ScanCoord%At1,ScanCoord%AT2, & |
| 790 | 1 | equemene | ScanCoord%At3,Geom,BprimT(1,I)) |
| 791 | 1 | equemene | CASE ('DIHEDRAL')
|
| 792 | 1 | equemene | CALL CONSTRAINTS_TORSION_DER2(Nat,ScanCoord%At1,ScanCoord%AT2, & |
| 793 | 1 | equemene | ScanCoord%At3,ScanCoord%At4,Geom,BprimT(1,I)) |
| 794 | 1 | equemene | ! CALL CONSTRAINTS_TORSION_DER(Nat,ScanCoord%At1,ScanCoord%AT2, & |
| 795 | 1 | equemene | ! ScanCoord%At3,ScanCoord%At4,Geom,BprimT(1,I)) |
| 796 | 1 | equemene | |
| 797 | 1 | equemene | END SELECT |
| 798 | 1 | equemene | ScanCoord => ScanCoord%next |
| 799 | 1 | equemene | END DO |
| 800 | 1 | equemene | |
| 801 | 1 | equemene | |
| 802 | 1 | equemene | !!!!!!!!!!!!!!!!!!!! |
| 803 | 1 | equemene | ! |
| 804 | 1 | equemene | ! PFL Debug |
| 805 | 1 | equemene | ! |
| 806 | 1 | equemene | !!!!!!!!!!!!!!!!!!!! |
| 807 | 1 | equemene | |
| 808 | 1 | equemene | if (debugPFL) THEN |
| 809 | 1 | equemene | WRITe(*,*) "DBG PFL Calc_Baker =====" |
| 810 | 1 | equemene | DzDc=0.d0 |
| 811 | 1 | equemene | WRITE(*,*) "PFL DBG, DzdC" |
| 812 | 1 | equemene | do I=1,Nat |
| 813 | 1 | equemene | Geom(1,i)=XyzGeomI(IGeomRef,1,orderInv(i)) ! XyzGeomI(NGeomI,3,Nat) |
| 814 | 1 | equemene | Geom(2,i)=XyzGeomI(IGeomRef,2,OrderInv(i)) |
| 815 | 1 | equemene | Geom(3,i)=XyzGeomI(IGeomRef,3,OrderInv(i)) |
| 816 | 1 | equemene | WRITE(*,*) I,OrderInv(I),Geom(:,I) |
| 817 | 1 | equemene | end do |
| 818 | 1 | equemene | CALL CalcBMat_int (nat, Geom, indzmat, dzdc, massat,atome) |
| 819 | 1 | equemene | |
| 820 | 1 | equemene | WRITE(*,*) "PFL DBG, BPrimT" |
| 821 | 1 | equemene | DO I=1,NPrim |
| 822 | 1 | equemene | WRITE(*,'(50(1X,F12.6))') BPrimT(:,I) |
| 823 | 1 | equemene | END DO |
| 824 | 1 | equemene | WRITe(*,*) "DBG PFL Calc_Baker =====" |
| 825 | 1 | equemene | END IF |
| 826 | 1 | equemene | |
| 827 | 1 | equemene | !!!!!!!!!!!!!!!!!!!! |
| 828 | 1 | equemene | ! |
| 829 | 1 | equemene | ! PFL Debug |
| 830 | 1 | equemene | ! |
| 831 | 1 | equemene | !!!!!!!!!!!!!!!!!!!! |
| 832 | 1 | equemene | |
| 833 | 1 | equemene | DEALLOCATE(Geom) |
| 834 | 1 | equemene | |
| 835 | 1 | equemene | ! BprimT(3*Nat,NPrim) |
| 836 | 1 | equemene | ! We now compute G=B(BT) matrix |
| 837 | 1 | equemene | ALLOCATE(Gmat(NPrim,NPrim)) |
| 838 | 1 | equemene | GMat=0.d0 |
| 839 | 1 | equemene | DO I=1,NPrim |
| 840 | 1 | equemene | DO J=1,3*Nat |
| 841 | 1 | equemene | GMat(:,I)=Gmat(:,I)+BprimT(J,:)*BprimT(J,I) !*1.d0/mass(atome(int(K/3.d0))) |
| 842 | 1 | equemene | END DO |
| 843 | 1 | equemene | END DO |
| 844 | 1 | equemene | |
| 845 | 1 | equemene | !!!!!!!!!!!!!!!!!!!! |
| 846 | 1 | equemene | ! |
| 847 | 1 | equemene | ! PFL Debug ---> |
| 848 | 1 | equemene | ! |
| 849 | 1 | equemene | !!!!!!!!!!!!!!!!!!!! |
| 850 | 1 | equemene | if (debugPFL) THEN |
| 851 | 1 | equemene | GMat=0.d0 |
| 852 | 1 | equemene | DO I=1,NPrim |
| 853 | 1 | equemene | GMat(I,I)=1. |
| 854 | 1 | equemene | END DO |
| 855 | 1 | equemene | END IF |
| 856 | 1 | equemene | |
| 857 | 1 | equemene | !!!!!!!!!!!!!!!!!!!! |
| 858 | 1 | equemene | ! |
| 859 | 1 | equemene | ! <--- PFL Debug |
| 860 | 1 | equemene | ! |
| 861 | 1 | equemene | !!!!!!!!!!!!!!!!!!!! |
| 862 | 1 | equemene | |
| 863 | 1 | equemene | !DO I=1,NPrim |
| 864 | 1 | equemene | ! WRITE(*,'(20(1X,F6.3))') GMat(1:min(20,NPrim),I) |
| 865 | 1 | equemene | !END DO |
| 866 | 1 | equemene | |
| 867 | 1 | equemene | ! We diagonalize G |
| 868 | 1 | equemene | ALLOCATE(EigVal(NPrim),EigVec(NPrim,NPrim)) |
| 869 | 1 | equemene | Call Jacobi(GMat,NPrim,EigVal,EigVec,NPrim) |
| 870 | 1 | equemene | Call Trie(NPrim,EigVal,EigVec,NPrim) |
| 871 | 1 | equemene | DO I=1,NPrim |
| 872 | 1 | equemene | WRITE(*,'(1X,"Vector ",I3,": e=",F8.3)') I,EigVal(i) |
| 873 | 1 | equemene | WRITE(*,'(20(1X,F8.4))') EigVec(1:min(20,NPrim),I) |
| 874 | 1 | equemene | END DO |
| 875 | 1 | equemene | |
| 876 | 1 | equemene | !DO I=1,NPrim |
| 877 | 1 | equemene | ! WRITE(*,'(1X,"Vector ",I3,(20(F8.3)))') I,EigVec(1:min(20,NPrim),I) |
| 878 | 1 | equemene | ! END DO |
| 879 | 1 | equemene | |
| 880 | 1 | equemene | !DO I=1,NPrim |
| 881 | 1 | equemene | ! WRITE(*,'(1X,"Vector ",I3,F8.3,F8.3,F8.3,F8.3,F8.3,F8.3,F8.3,F8.3,F8.3,F8.3,F8.3,F8.3,& |
| 882 | 1 | equemene | ! F8.3,F8.3,F8.3,F8.3,F8.3)') I,EigVec(1,I),EigVec(4,I),EigVec(6,I),EigVec(14,I),EigVec(16,I),& |
| 883 | 1 | equemene | !EigVec(2,I),EigVec(3,I),EigVec(5,I),EigVec(12,I),EigVec(13,I),EigVec(15,I),EigVec(7,I),& |
| 884 | 1 | equemene | !EigVec(8,I),EigVec(9,I),EigVec(10,I),EigVec(11,I),EigVec(17,I) |
| 885 | 1 | equemene | !END DO |
| 886 | 1 | equemene | |
| 887 | 1 | equemene | ! We construct the new DzDc(b=baker) matrix |
| 888 | 1 | equemene | ! UMat is nonredundant vector set, i.e. set of eigenvectors of BB^T corresponding |
| 889 | 1 | equemene | ! to eigenvalues > zero. |
| 890 | 1 | equemene | ALLOCATE(UMat(NPrim,NCoord)) |
| 891 | 1 | equemene | ALLOCATE(UMat_local(NPrim,NCoord)) |
| 892 | 1 | equemene | ALLOCATE(UMat_tmp(NPrim,3*Nat-6)) |
| 893 | 1 | equemene | ! BMat_BakerT(3*Nat,NCoord), allocated in Path.f90, |
| 894 | 1 | equemene | ! NCoord=3*Nat-6 |
| 895 | 1 | equemene | BMat_BakerT = 0.d0 |
| 896 | 1 | equemene | J=0 |
| 897 | 1 | equemene | UMat=0.d0 |
| 898 | 1 | equemene | UMat_tmp=0.d0 |
| 899 | 1 | equemene | DO I=1,NPrim |
| 900 | 1 | equemene | IF (abs(eigval(I))>=1e-6) THEN |
| 901 | 1 | equemene | J=J+1 |
| 902 | 1 | equemene | DO K=1,NPrim |
| 903 | 1 | equemene | !DzDcb(:,J)=DzDcb(:,J)+Eigvec(K,I)*BprimT(:,K) ! original |
| 904 | 1 | equemene | ! BprimT is transpose of B^prim. |
| 905 | 1 | equemene | ! B = UMat^T B^prim, B^T = (B^prim)^T UMat |
| 906 | 1 | equemene | BMat_BakerT(:,J)=BMat_BakerT(:,J)+BprimT(:,K)*Eigvec(K,I) |
| 907 | 1 | equemene | !Dzdcb(:,J)=DzDcb(:,J)+Eigvec(K,:)*BprimT(I,K) |
| 908 | 1 | equemene | END DO |
| 909 | 1 | equemene | IF(J .GT. 3*Nat-6) THEN |
| 910 | 1 | equemene | WRITE(*,*) 'Number of vectors in Eigvec with eigval .GT. 1e-6(=UMat) (=' & |
| 911 | 1 | equemene | ,J,') exceeded 3*Nat-6=',3*Nat-6, & |
| 912 | 1 | equemene | 'Stopping the calculation.' |
| 913 | 1 | equemene | STOP |
| 914 | 1 | equemene | END IF |
| 915 | 1 | equemene | UMat_tmp(:,J) = Eigvec(:,I) |
| 916 | 1 | equemene | END IF |
| 917 | 1 | equemene | END DO |
| 918 | 1 | equemene | |
| 919 | 1 | equemene | if (debug) THEN |
| 920 | 1 | equemene | WRITE(*,*) "DBG Calc_Baker: UMat" |
| 921 | 1 | equemene | DO I=1,3*Nat-6 |
| 922 | 1 | equemene | WRITE(*,'(50(1X,F12.8))') UMat_tmp(:,I) |
| 923 | 1 | equemene | END DO |
| 924 | 1 | equemene | END IF |
| 925 | 1 | equemene | |
| 926 | 1 | equemene | ! THIS IS TO BE CORRECTED: |
| 927 | 1 | equemene | UMat = UMat_tmp |
| 928 | 1 | equemene | !J=0 |
| 929 | 1 | equemene | !DO I=1, 3*Nat-6 |
| 930 | 1 | equemene | ! IF ((I .NE. 6) .AND. (I .NE. 7) .AND. (I .NE. 9)) Then |
| 931 | 1 | equemene | ! J=J+1 |
| 932 | 1 | equemene | ! Print *, 'J=', J, ' I=', I |
| 933 | 1 | equemene | ! UMat(:,J) = UMat_tmp(:,I) |
| 934 | 1 | equemene | !END IF |
| 935 | 1 | equemene | !END DO |
| 936 | 1 | equemene | ! BMat_BakerT=0.d0 |
| 937 | 1 | equemene | !DO I=1,NCoord |
| 938 | 1 | equemene | ! DO J=1,NPrim |
| 939 | 1 | equemene | ! B = UMat^T B^prim, B^T = (B^prim)^T UMat |
| 940 | 1 | equemene | ! BMat_BakerT(:,I)=BMat_BakerT(:,I)+BprimT(:,J)*UMat(J,I) |
| 941 | 1 | equemene | ! END DO |
| 942 | 1 | equemene | !END DO |
| 943 | 1 | equemene | ! TO BE CORRECTED, ENDS HERE. |
| 944 | 1 | equemene | |
| 945 | 1 | equemene | IntCoordI=0.d0 |
| 946 | 1 | equemene | DO IGeom=1, NGeomI |
| 947 | 1 | equemene | DO I=1, NPrim |
| 948 | 1 | equemene | ! Transpose of UMat is needed below, that is why UMat(I,:). |
| 949 | 1 | equemene | IntCoordI(IGeom,:) = IntCoordI(IGeom,:) + UMat(I,:)*Xprimitive(IGeom,I) |
| 950 | 1 | equemene | END DO |
| 951 | 1 | equemene | END DO |
| 952 | 1 | equemene | |
| 953 | 1 | equemene | if (debug) THEN |
| 954 | 1 | equemene | WRITE(*,*) "DBG Calc_Baker IntCoordI" |
| 955 | 1 | equemene | DO IGeom=1, NGeomI |
| 956 | 1 | equemene | WRITE(*,'(1X,I5,":",30(1X,F10.6))') IGeom,IntCoordI(IGeom,:) |
| 957 | 1 | equemene | END DO |
| 958 | 1 | equemene | END IF |
| 959 | 1 | equemene | |
| 960 | 1 | equemene | |
| 961 | 1 | equemene | ! the following might be used to get rid of some coordinates that |
| 962 | 1 | equemene | ! do not respect the symmetry of the problem. |
| 963 | 1 | equemene | |
| 964 | 1 | equemene | ! ! We look at the group symmetry of the molecule. Nothing complicate here |
| 965 | 1 | equemene | ! ! we just look for planar symmetry. |
| 966 | 1 | equemene | ! ! We first place it into the standard orientation |
| 967 | 1 | equemene | ! Call Std_ori(Nat,x,y,z,a,b,c,massat) |
| 968 | 1 | equemene | |
| 969 | 1 | equemene | ! ! Is the molecule planar |
| 970 | 1 | equemene | ! DO I=1,Nat |
| 971 | 1 | equemene | ! WRITE(*,*) nom(atome(i)),x(i),y(i),z(i) |
| 972 | 1 | equemene | ! END DO |
| 973 | 1 | equemene | DEALLOCATE(BprimT,GMat,EigVal,EigVec) ! BprimT is B^prim^T |
| 974 | 1 | equemene | ! Calculation of BTransInv starts here: |
| 975 | 1 | equemene | ! Calculation of BBT(3*Nat-6,3*Nat-6)=BB^T: |
| 976 | 1 | equemene | ! BMat_BakerT(3*Nat,NCoord) is Transpose of B = UMat^TB^prim |
| 977 | 1 | equemene | ALLOCATE(BBT(NCoord,NCoord)) |
| 978 | 1 | equemene | BBT = 0.d0 |
| 979 | 1 | equemene | DO I=1, NCoord |
| 980 | 1 | equemene | DO J=1, 3*Nat |
| 981 | 1 | equemene | ! BBT(:,I) forms BB^T |
| 982 | 1 | equemene | BBT(:,I) = BBT(:,I) + BMat_BakerT(J,:)*BMat_BakerT(J,I) |
| 983 | 1 | equemene | END DO |
| 984 | 1 | equemene | END DO |
| 985 | 1 | equemene | |
| 986 | 1 | equemene | !ALLOCATE(V(3*Nat,3*Nat)) |
| 987 | 1 | equemene | !Call GINVSE(BBT,3*Nat,3*Nat,BBT_inv,3*Nat,3*Nat,3*Nat,1e-6,V,3*Nat,FAIL,NOUT) |
| 988 | 1 | equemene | !DEALLOCATE(V) |
| 989 | 1 | equemene | ALLOCATE(BBT_inv(NCoord,NCoord)) |
| 990 | 1 | equemene | ! GenInv is in Mat_util.f90 |
| 991 | 1 | equemene | Call GenInv(NCoord,BBT,BBT_inv,NCoord) |
| 992 | 1 | equemene | |
| 993 | 1 | equemene | ! Calculation of (B^T)^-1 = (BB^T)^-1B: |
| 994 | 1 | equemene | !ALLOCATE(BTransInv(3*Nat-6,3*Nat)) ! allocated in Path.f90 |
| 995 | 1 | equemene | BTransInv = 0.d0 |
| 996 | 1 | equemene | DO I=1, 3*Nat |
| 997 | 1 | equemene | DO J=1, NCoord |
| 998 | 1 | equemene | BTransInv(:,I) = BTransInv(:,I) + BBT_inv(:,J)*BMat_BakerT(I,J) |
| 999 | 1 | equemene | END DO |
| 1000 | 1 | equemene | END DO |
| 1001 | 1 | equemene | BTransInv_local = BTransInv |
| 1002 | 1 | equemene | UMat_local = UMat |
| 1003 | 1 | equemene | DEALLOCATE(BBT,BBT_inv,UMat_tmp) |
| 1004 | 1 | equemene | DEALLOCATE(x_refGeom,y_refGeom,z_refGeom,x,y,z) |
| 1005 | 1 | equemene | IF (debug) WRITE(*,*) "DBG Calc_baker over." |
| 1006 | 1 | equemene | END SUBROUTINE Calc_baker |