Statistiques
| Révision :

root / src / Calc_Xprim.f90 @ 7

Historique | Voir | Annoter | Télécharger (6,74 ko)

1
SUBROUTINE Calc_Xprim(nat,x,y,z,Coordinate,NPrim,XPrimitive,XPrimRef)
2
!     
3
! This subroutine uses the description of a list of Coordinates
4
! to compute the values of the coordinates for a given geometry.
5
!
6
!!!!!!!!!!
7
! Input:
8
!   Na: INTEGER, Number of atoms
9
!  x,y,z(Na): REAL, cartesian coordinates of the considered geometry
10
! Coordinate (Pointer(ListCoord)): description of the wanted coordiantes
11
! NPrim, INTEGER: Number of coordinates to compute
12
!
13
! Optional: XPrimRef(NPrim) REAL: array that contains coordinates values for
14
! a former geometry. Useful for Dihedral angles evolution...
15

    
16
!!!!!!!!!!!
17
! Output:
18
! XPrimimite(NPrim) REAL: array that will contain the values of the coordinates
19
!
20
!!!!!!!!!
21

    
22
  Use VarTypes
23
  Use Io_module
24
  Use Path_module, only : pi
25

    
26
  IMPLICIT NONE
27

    
28
  Type (ListCoord), POINTER :: Coordinate
29
  INTEGER(KINT), INTENT(IN) :: Nat,NPrim
30
  REAL(KREAL), INTENT(IN) :: x(Nat), y(Nat), z(Nat)
31
  REAL(KREAL), INTENT(IN), OPTIONAL :: XPrimRef(NPrim) 
32
  REAL(KREAL), INTENT(OUT) :: XPrimitive(NPrim)
33

    
34

    
35
  Type (ListCoord), POINTER :: ScanCoord
36

    
37
  real(KREAL) :: vx1,vy1,vz1,norm1
38
  real(KREAL) :: vx2,vy2,vz2,norm2
39
  real(KREAL) :: vx3,vy3,vz3,norm3
40
  real(KREAL) :: vx4,vy4,vz4,norm4
41
  real(KREAL) :: vx5,vy5,vz5,norm5
42

    
43
  INTEGER(KINT) :: I
44

    
45
  REAL(KREAL) :: DiheTmp
46

    
47
  LOGICAL :: debug, debugPFL
48

    
49
  INTERFACE
50
     function valid(string) result (isValid)
51
       CHARACTER(*), intent(in) :: string
52
       logical                  :: isValid
53
     END function VALID
54

    
55
     FUNCTION angle(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2)
56
        use Path_module, only :  Pi,KINT, KREAL
57
       real(KREAL) ::  v1x,v1y,v1z,norm1
58
       real(KREAL) ::  v2x,v2y,v2z,norm2
59
       real(KREAL) ::  angle
60
     END FUNCTION ANGLE
61

    
62
     FUNCTION SinAngle(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2)
63
        use Path_module, only :  Pi,KINT, KREAL
64
       real(KREAL) ::  v1x,v1y,v1z,norm1
65
       real(KREAL) ::  v2x,v2y,v2z,norm2
66
       real(KREAL) ::  SinAngle
67
     END FUNCTION SINANGLE
68

    
69

    
70
     FUNCTION angle_d(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2,v3x,v3y,v3z,norm3)
71
        use Path_module, only :  Pi,KINT, KREAL
72
       real(KREAL) ::  v1x,v1y,v1z,norm1
73
       real(KREAL) ::  v2x,v2y,v2z,norm2
74
       real(KREAL) ::  v3x,v3y,v3z,norm3
75
       real(KREAL) ::  angle_d,ca,sa
76
     END FUNCTION ANGLE_D
77

    
78
  END INTERFACE
79

    
80

    
81
  debug=valid("Calc_Xprim")
82
  debugPFL=valid("BakerPFL")
83

    
84
  if (debug) Call Header("Entering Cal_XPrim")
85

    
86

    
87
  IF (debug) THEN
88
     WRITE(*,*) "DBG Calc_Xprim, geom"
89
     DO I=1,Nat
90
        WRITE(*,'(1X,I5,3(1X,F15.3))') I,X(I),y(i),z(i)
91
     END DO
92
     IF (Present(XPrimRef)) THEN
93
        WRITE(*,*) "XPrimRef"
94
        WRITE(*,'(15(1X,F10.6))') XPrimRef
95
     END IF
96
     WRITE(*,*) "NPrim:",NPrim
97
  END IF
98

    
99

    
100
  WRITE(*,*) "Coordinate:",associated(Coordinate),coordinate%Type
101

    
102
     ScanCoord => Coordinate
103

    
104
  WRITE(*,*) "ScanCoord:",associated(ScanCoord),ScanCoord%Type
105

    
106
!     WRITE(*,*) "coucou"
107
     I=0 ! index for the NPrim (NPrim is the number of primitive coordinates).
108
     DO WHILE (Associated(ScanCoord%next))
109
        I=I+1
110
!        WRITE(*,*) i
111
        SELECT CASE (ScanCoord%Type)
112
        CASE ('BOND')
113
           Call vecteur(ScanCoord%At2,ScanCoord%At1,x,y,z,vx2,vy2,vz2,Norm2)
114
           Xprimitive(I) = Norm2  
115
        CASE ('ANGLE')
116
           Call vecteur(ScanCoord%At2,ScanCoord%At3,x,y,z,vx1,vy1,vz1,Norm1)
117
           Call vecteur(ScanCoord%At2,ScanCoord%At1,x,y,z,vx2,vy2,vz2,Norm2)
118
           Xprimitive(I) = angle(vx1,vy1,vz1,Norm1,vx2,vy2,vz2,Norm2)*Pi/180.
119
        CASE ('DIHEDRAL')
120
           Call vecteur(ScanCoord%At2,ScanCoord%At1,x,y,z,vx1,vy1,vz1,Norm1)
121
           Call vecteur(ScanCoord%At2,ScanCoord%At3,x,y,z,vx2,vy2,vz2,Norm2)
122
           Call vecteur(ScanCoord%At3,ScanCoord%At4,x,y,z,vx3,vy3,vz3,Norm3)
123
           Call produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, &
124
                vx4,vy4,vz4,norm4)
125
           Call produit_vect(vx3,vy3,vz3,vx2,vy2,vz2, &
126
                vx5,vy5,vz5,norm5)
127
           
128
           DiheTmp= angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, &
129
                vx2,vy2,vz2,norm2)
130
           Xprimitive(I) = DiheTmp*Pi/180.
131
! PFL 15th March 2008 ->
132
! I think that the test on changes less than Pi should be enough
133
!! We treat large dihedral angles differently as +180=-180 mathematically and physically
134
!! but this causes lots of troubles when converting baker to cart.
135
!! So we ensure that large dihedral angles always have the same sign
136
!                    if (abs(ScanCoord%SignDihedral).NE.1) THEN
137
!                        ScanCoord%SignDihedral=Int(Sign(1.D0,DiheTmp))
138
!                    ELSE
139
!                        If ((abs(DiheTmp).GE.170.D0).AND.(Sign(1.,DiheTmp)*ScanCoord%SignDihedral<0)) THEN
140
!                          Xprimitive(I) = DiheTmp*Pi/180.+ ScanCoord%SignDihedral*2.*Pi
141
!                       END IF
142
!                  END IF
143
!!!! <- PFL 15th March 2008
144
! Another test... will all this be consistent ??? 
145
! We want the shortest path, so we check that the change in dihedral angles is less than Pi:
146
           IF (Present(XPrimRef)) THEN
147
              IF (abs(Xprimitive(I)-XPrimRef(I)).GE.Pi) THEN
148
                 if (debug) THEN 
149
                    WRITE(*,*) "Pb dihedral, i,Xprimivite,XPrimref=",i,XPrimitive(I),XPrimRef(I)
150
                    WRITE(*,*) "In deg Xprimivite,XPrimref=",XPrimitive(I)*180./Pi,XPrimRef(I)*180/Pi
151
                 END IF
152
                 if ((Xprimitive(I)-XPrimRef(I)).GE.Pi) THEN
153
                    Xprimitive(I)=Xprimitive(I)-2.d0*Pi
154
                 ELSE
155
                    Xprimitive(I)=Xprimitive(I)+2.d0*Pi
156
                 END IF
157
              END IF
158
              if (debug) WRITE(*,*) " New Xprimivite",XPrimitive(I),XPrimitive(I)*180./Pi
159
           END IF
160
        CASE DEFAULT
161
           WRITE(*,*) "Scancoord%type unknown for I",I,scancoord%type
162
      END SELECT
163
        ScanCoord => ScanCoord%next
164
     END DO ! matches DO WHILE (Associated(ScanCoord%next))
165

    
166

    
167

    
168
  IF (debug) THEN
169
     WRITE(*,*) "DBG Calc_Xprim Values"
170
     WRITE(*,*) "Found ",I," primitives"
171
     
172
     ScanCoord=>Coordinate
173
     I=0 ! index for the NPrim (NPrim is the number of primitive coordinates).
174
     DO WHILE (Associated(ScanCoord%next))
175
        I=I+1
176
        SELECT CASE (ScanCoord%Type)
177
        CASE ('BOND')
178
           WRITE(*,'(1X,I3,":",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1,ScanCoord%At2,Xprimitive(I)
179
        CASE ('ANGLE')
180
           WRITE(*,'(1X,I3,":",I5," - ",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1, &
181
            ScanCoord%At2, ScanCoord%At3,Xprimitive(I)*180./Pi
182
        CASE ('DIHEDRAL')
183
           WRITE(*,'(1X,I3,":",I5," - ",I5," - ",I5," - ",I5," = ",F15.6)') I,ScanCoord%At1,&
184
           ScanCoord%At2, ScanCoord%At3,ScanCoord%At4,Xprimitive(I)*180./Pi
185
        END SELECT
186
        ScanCoord => ScanCoord%next
187
     END DO ! matches DO WHILE (Associated(ScanCoord%next))
188
  END IF
189

    
190
  if (debug) Call Header(" Cal_XPrim OVER ")
191

    
192
END SUBROUTINE Calc_Xprim