Révision 7 src/Calc_Xprim.f90

Calc_Xprim.f90 (revision 7)
80 80

  
81 81
  debug=valid("Calc_Xprim")
82 82
  debugPFL=valid("BakerPFL")
83
  if (debug) WRITE(*,*) "============= Entering Cal_XPrim =============="
84 83

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

  
86

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

  
95 99

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

  
98
     ScanCoord=>Coordinate
102
     ScanCoord => Coordinate
99 103

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

  
102 106
!     WRITE(*,*) "coucou"
103 107
     I=0 ! index for the NPrim (NPrim is the number of primitive coordinates).
......
120 124
                vx4,vy4,vz4,norm4)
121 125
           Call produit_vect(vx3,vy3,vz3,vx2,vy2,vz2, &
122 126
                vx5,vy5,vz5,norm5)
123

  
124
                 DiheTmp= angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, &
125
                         vx2,vy2,vz2,norm2)
126
                 Xprimitive(I) = DiheTmp*Pi/180.
127
           
128
           DiheTmp= angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, &
129
                vx2,vy2,vz2,norm2)
130
           Xprimitive(I) = DiheTmp*Pi/180.
127 131
! PFL 15th March 2008 ->
128 132
! I think that the test on changes less than Pi should be enough
129 133
!! We treat large dihedral angles differently as +180=-180 mathematically and physically
......
139 143
!!!! <- PFL 15th March 2008
140 144
! Another test... will all this be consistent ??? 
141 145
! We want the shortest path, so we check that the change in dihedral angles is less than Pi:
142
                  IF (Present(XPrimRef)) THEN
143
                     IF (abs(Xprimitive(I)-XPrimRef(I)).GE.Pi) THEN
144
                       if (debug) THEN 
145
                         WRITE(*,*) "Pb dihedral, i,Xprimivite,XPrimref=",i,XPrimitive(I),XPrimRef(I)
146
                         WRITE(*,*) "In deg Xprimivite,XPrimref=",XPrimitive(I)*180./Pi,XPrimRef(I)*180/Pi
147
                       END IF
148
                        if ((Xprimitive(I)-XPrimRef(I)).GE.Pi) THEN
149
                           Xprimitive(I)=Xprimitive(I)-2.d0*Pi
150
                        ELSE
151
                           Xprimitive(I)=Xprimitive(I)+2.d0*Pi
152
                        END IF
153
                     END IF
154
                       if (debug) WRITE(*,*) " New Xprimivite",XPrimitive(I),XPrimitive(I)*180./Pi
155
                  END IF
156
        END SELECT
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
157 163
        ScanCoord => ScanCoord%next
158 164
     END DO ! matches DO WHILE (Associated(ScanCoord%next))
159 165

  
......
161 167

  
162 168
  IF (debug) THEN
163 169
     WRITE(*,*) "DBG Calc_Xprim Values"
170
     WRITE(*,*) "Found ",I," primitives"
164 171
     
165 172
     ScanCoord=>Coordinate
166 173
     I=0 ! index for the NPrim (NPrim is the number of primitive coordinates).
......
179 186
        ScanCoord => ScanCoord%next
180 187
     END DO ! matches DO WHILE (Associated(ScanCoord%next))
181 188
  END IF
182
  if (debug) WRITE(*,*) "============= Cal_XPrim OVER =============="
189

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

  
183 192
END SUBROUTINE Calc_Xprim

Formats disponibles : Unified diff