Statistiques
| Révision :

root / src / Calc_baker_allGeomF.f90 @ 8

Historique | Voir | Annoter | Télécharger (7,31 ko)

1
SUBROUTINE Calc_baker_allGeomF()
2
  !     
3
  ! This subroutine analyses a geometry to construct the baker
4
  ! delocalized internal coordinates
5
  ! v1.0
6
  ! We use only one geometry
7
  !     
8
  Use Path_module, only : BMat_BakerT,Nat,NCoord,UMatF, &
9
                          NPrim,BTransInvF,Coordinate, &
10
              ScanCoord,BprimT,BBT,BBT_inv,XprimitiveF, &
11
              NgeomF,XyzGeomF
12
  ! BMat_BakerT(3*Nat,NCoord), NCoord=3*Nat or NFree=3*Nat-6-Symmetry_elimination
13
  ! depending upon the coordinate choice. IntCoordI(NGeomI,NCoord) where 
14
  ! UMatF(NGeomI,NPrim,NCoord), NCoord number of vectors in UMat matrix, i.e. NCoord
15
  ! Baker coordinates. NPrim is the number of primitive internal coordinates.
16
  
17
  Use Io_module
18
  IMPLICIT NONE
19

    
20
  REAL(KREAL), ALLOCATABLE :: Geom(:,:) !(3,Nat)
21
  ! NPrim is the number of primitive coordinates and NCoord is the number
22
  ! of internal coordinates. BMat is actually (NPrim,3*Nat).
23
  REAL(KREAL), ALLOCATABLE :: GMat(:,:) !(NPrim,NPrim)
24
  ! EigVec(..) contains ALL eigevectors of BMat times BprimT, NOT only Baker Coordinate vectors.
25
  REAL(KREAL), ALLOCATABLE :: EigVec(:,:), EigVal(:) ! EigVec(NPrim,NPrim)
26
  REAL(KREAL), ALLOCATABLE :: x(:), y(:), z(:)
27
  REAL(KREAL), ALLOCATABLE :: XPrimRef(:) ! NPrim
28
  INTEGER(KINT) :: IGeom
29
  
30

    
31
  INTEGER(KINT) :: I, J, K
32

    
33

    
34
  LOGICAL :: debug
35
  LOGICAL :: DebugPFL
36

    
37
  INTERFACE
38
     function valid(string) result (isValid)
39
       CHARACTER(*), intent(in) :: string
40
       logical                  :: isValid
41
     END function VALID
42

    
43
     FUNCTION angle(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2)
44
       use Path_module, only :  Pi,KINT, KREAL
45
       real(KREAL) ::  v1x,v1y,v1z,norm1
46
       real(KREAL) ::  v2x,v2y,v2z,norm2
47
       real(KREAL) ::  angle
48
     END FUNCTION ANGLE
49

    
50
     FUNCTION angle_d(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2,v3x,v3y,v3z,norm3)
51
       use Path_module, only :  Pi,KINT, KREAL
52
       real(KREAL) ::  v1x,v1y,v1z,norm1
53
       real(KREAL) ::  v2x,v2y,v2z,norm2
54
       real(KREAL) ::  v3x,v3y,v3z,norm3
55
       real(KREAL) ::  angle_d,ca,sa
56
     END FUNCTION ANGLE_D
57

    
58

    
59

    
60
     SUBROUTINE Calc_Xprim(nat,x,y,z,Coordinate,NPrim,XPrimitive,XPrimRef)
61
       !     
62
       ! This subroutine uses the description of a list of Coordinates
63
       ! to compute the values of the coordinates for a given geometry.
64
       !
65
!!!!!!!!!!
66
       ! Input:
67
       !   Na: INTEGER, Number of atoms
68
       !  x,y,z(Na): REAL, cartesian coordinates of the considered geometry
69
       ! Coordinate (Pointer(ListCoord)): description of the wanted coordiantes
70
       ! NPrim, INTEGER: Number of coordinates to compute
71
       !
72
       ! Optional: XPrimRef(NPrim) REAL: array that contains coordinates values for
73
       ! a former geometry. Useful for Dihedral angles evolution...
74

    
75
!!!!!!!!!!!
76
       ! Output:
77
       ! XPrimimite(NPrim) REAL: array that will contain the values of the coordinates
78
       !
79
!!!!!!!!!
80

    
81
       Use VarTypes
82
       Use Io_module
83
       Use Path_module, only : pi
84

    
85
       IMPLICIT NONE
86

    
87
       Type (ListCoord), POINTER :: Coordinate
88
       INTEGER(KINT), INTENT(IN) :: Nat,NPrim
89
       REAL(KREAL), INTENT(IN) :: x(Nat), y(Nat), z(Nat)
90
       REAL(KREAL), INTENT(IN), OPTIONAL :: XPrimRef(NPrim) 
91
       REAL(KREAL), INTENT(OUT) :: XPrimitive(NPrim)
92

    
93
     END SUBROUTINE CALC_XPRIM
94
  END INTERFACE
95

    
96

    
97

    
98

    
99
  debug=valid("Calc_baker_allGeomF")
100
  debugPFL=valid("bakerPFL")
101
  if (debug) WRITE(*,*) '============ Entering Calc_baker_allGeomF ============='
102

    
103
  ALLOCATE(Geom(3,Nat),x(Nat),y(Nat),z(Nat))
104
  ALLOCATE(XPrimRef(NPrim))
105
 
106
  ! Now calculating values of all primitive bonds for all final geometries:
107
  DO IGeom=1, NGeomF
108
     x(1:Nat) = XyzGeomF(IGeom,1,1:Nat)
109
     y(1:Nat) = XyzGeomF(IGeom,2,1:Nat)
110
     z(1:Nat) = XyzGeomF(IGeom,3,1:Nat)
111
     XPrimREf=XPrimitiveF(IGeom,:)
112
     Call Calc_XPrim(nat,x,y,z,Coordinate,NPrim,XPrimitiveF(IGeom,:),XPrimRef)
113
  END DO ! matches DO IGeom=1, NGeomF
114
   
115
  ALLOCATE(BprimT(3*Nat,NPrim))
116
  ALLOCATE(Gmat(NPrim,NPrim))
117
  ALLOCATE(EigVal(NPrim),EigVec(NPrim,NPrim))
118
  ALLOCATE(BBT(NCoord,NCoord))
119
  ALLOCATE(BBT_inv(NCoord,NCoord))
120
  BTransInvF = 0.d0
121

    
122
  DO IGeom=1, NGeomF 
123
     Geom(1,:)=XyzGeomF(IGeom,1,1:Nat) ! XyzGeomI(NGeomI,3,Nat)
124
     Geom(2,:)=XyzGeomF(IGeom,2,1:Nat)
125
     Geom(3,:)=XyzGeomF(IGeom,3,1:Nat)
126

    
127
   BprimT=0.d0
128
     ScanCoord=>Coordinate
129
     I=0
130
     DO WHILE (Associated(ScanCoord%next))
131
        I=I+1
132
        SELECT CASE (ScanCoord%Type)
133
        CASE ('BOND') 
134
          CALL CONSTRAINTS_BONDLENGTH_DER(Nat,ScanCoord%at1,ScanCoord%AT2, &
135
               Geom,BprimT(1,I))
136
        CASE ('ANGLE')
137
          CALL CONSTRAINTS_BONDANGLE_DER(Nat,ScanCoord%At1,ScanCoord%AT2, &
138
               ScanCoord%At3,Geom,BprimT(1,I))
139
        CASE ('DIHEDRAL')
140
          CALL CONSTRAINTS_TORSION_DER2(Nat,ScanCoord%At1,ScanCoord%AT2, &
141
               ScanCoord%At3,ScanCoord%At4,Geom,BprimT(1,I))
142
      END SELECT
143
        ScanCoord => ScanCoord%next
144
     END DO
145

    
146
   ! BprimT(3*Nat,NPrim)
147
   ! We now compute G=B(BT) matrix
148
   GMat=0.d0
149
   DO I=1,NPrim
150
    DO J=1,3*Nat
151
           GMat(:,I)=Gmat(:,I)+BprimT(J,:)*BprimT(J,I) !*1.d0/mass(atome(int(K/3.d0)))
152
    END DO
153
   END DO
154
  
155
   ! Diagonalize G
156
   EigVal=0.d0
157
   EigVec=0.d0
158
   Call Jacobi(GMat,NPrim,EigVal,EigVec,NPrim)
159
   Call Trie(NPrim,EigVal,EigVec,NPrim)
160
   DO I=1,NPrim
161
    !WRITE(*,'(1X,"Vector ",I3,": e=",F8.3)') I,EigVal(i)
162
    !WRITE(*,'(20(1X,F8.4))') EigVec(1:min(20,NPrim),I)
163
   END DO
164

    
165
   ! UMatF is nonredundant vector set, i.e. set of eigenvectors of BB^T
166
   !  corresponding to eigenvalues > zero.
167
   ! BMat_BakerT(3*Nat,NCoord), allocated in Path.f90, 
168
   ! NCoord=3*Nat-6
169
   BMat_BakerT = 0.d0
170
   J=0
171
   DO I=1,NPrim
172
    IF (abs(eigval(I))>=1e-6) THEN
173
       J=J+1
174
       DO K=1,NPrim
175
        ! BprimT is transpose of B^prim.
176
        ! B = UMatF^T B^prim, B^T = (B^prim)^T UMatF
177
        BMat_BakerT(:,J)=BMat_BakerT(:,J)+BprimT(:,K)*Eigvec(K,I)
178
       END DO
179
         IF(J .GT. 3*Nat-6) THEN
180
           WRITE(*,*) 'Number of vectors in Eigvec with eigval .GT. 1e-6(=UMatF) (=' &
181
                      ,J,') exceeded 3*Nat-6=',3*Nat-6, &
182
            'Stopping the calculation.'
183
           STOP
184
         END IF
185
         UMatF(IGeom,:,J) =  Eigvec(:,I)
186
    END IF
187
   END DO
188

    
189
!!!!!!!!!!!!!!!!!!!!
190
!
191
! Debug purposes
192
!
193
         if (debugPFL) THEN
194
        UMatF(IGeom,:,:)=0.
195
        DO J=1,3*Nat-6
196
           UMatF(IGeom,J,J)=1.
197
        END DO
198
        END IF
199

    
200
 
201
   !DO I=1, NPrim ! This loop is not needed because we already have IntCoordF
202
             ! from interpolation. 
203
   ! Transpose of UMatF is needed below, that is why UMatF(IGeom,I,:).
204
   ! IntCoordF(IGeom,:) = IntCoordF(IGeom,:) + UMat(IGeom,I,:)*XprimitiveF(IGeom,I)
205
   !END DO
206

    
207
   ! Calculation of BTransInvF starts here:
208
   ! Calculation of BBT(3*Nat-6,3*Nat-6)=BB^T:
209
   ! BMat_BakerT(3*Nat,NCoord) is Transpose of B = UMatF^TB^prim
210
    
211
   BBT = 0.d0
212
   DO I=1, NCoord
213
      DO J=1, 3*Nat
214
       ! BBT(:,I) forms BB^T
215
       BBT(:,I) = BBT(:,I) + BMat_BakerT(J,:)*BMat_BakerT(J,I)
216
    END DO
217
   END DO
218
    
219
   Call GenInv(NCoord,BBT,BBT_inv,NCoord) ! GenInv is in Mat_util.f90
220

    
221
   ! Calculation of (B^T)^-1 = (BB^T)^-1B:
222
   DO I=1, 3*Nat
223
      DO J=1, NCoord
224
        BTransInvF(IGeom,:,I) = BTransInvF(IGeom,:,I) + BBT_inv(:,J)*BMat_BakerT(I,J)
225
      END DO
226
   END DO
227
   
228
  END DO !matches DO IGeom=1, NGeomF
229
  
230
  DEALLOCATE(BBT,BBT_inv,BprimT,GMat,EigVal,EigVec)
231
  DEALLOCATE(Geom,x,y,z,XprimRef)
232
  
233
  IF (debug) WRITE(*,*) "DBG Calc_baker_allGeomF over."
234
END SUBROUTINE Calc_baker_allGeomF