root / src / Egrad.f90 @ 4
Historique  Voir  Annoter  Télécharger (15,19 ko)
1 
subroutine Egrad(E,Geom,Grad,NCoord,IGeom,IOpt,GeomCart,FOptGeom,GeomOld,GeomCart_old) 

2  
3 
! This routines calculates the energy E and the gradient Grad of 
4 
! a molecule with Geometry Geom (may be in internal coordinates), 
5 
! using for now, either Gaussian or Ext, more general later. 
6  
7 
use Path_module, only : Nat,AtName,Coord,dzdc,indzmat,Nom,Atome,massat,unit, & 
8 
prog,NCart,XyzGeomF,IntCoordF,BTransInv, XyzGeomI, renum, & 
9 
GeomOld_all,BTransInvF,BTransInv_local,UMatF,UMat_local & 
10 
, BprimT,a0,ScanCoord, Coordinate,NPrim,BMat_BakerT,BBT,BBT_inv & 
11 
, Order,OrderInv, XPrimitiveF 
12  
13 
! IntCoordF(NGeomF,NCoord),GeomOld_all(NGeomF,NCoord) 
14 
! allocated in Path.f90 
15  
16 
use Io_module 
17 
IMPLICIT NONE 
18  
19 
! Energy (calculated if F300K=.F., else estimated) 
20 
REAL(KREAL), INTENT (OUT) :: E 
21 
! NCoord: Number of the degrees of freedom 
22 
! IGeom: index of the geometry. 
23 
INTEGER(KINT), INTENT (IN) :: NCoord, IGeom, IOpt 
24 
! Geometry at which gradient is calculated (cf Factual also): 
25 
REAL(KREAL), INTENT (INOUT) :: Geom(NCoord) 
26 
! Gradient calculated at Geom geometry: 
27 
REAL(KREAL), INTENT (OUT) :: Grad(NCoord) 
28 
! Cartesian geometry corresponding to (Internal Geometry) Geom: 
29 
REAL(KREAL), INTENT (OUT) :: GeomCart(Nat,3) 
30 
!!! Optional, just for geometry optimization with Baker coordinates 
31 
REAL(KREAL), INTENT (IN), OPTIONAL :: GeomCart_old(Nat,3) 
32 
REAL(KREAL), INTENT (INOUT), OPTIONAL :: GeomOld(NCoord) 
33 
! FOptGeom is a flag indicating if we are doing a geom optimization 
34 
! it can be omitted so that we use a local flag: Flag_Opt_Geom 
35 
LOGICAL,INTENT (IN), OPTIONAL :: FOptGeom 
36 
! if FOptGeom is given Flag_Opt_Geom=FOptGeom, else Flag_Opt_Geom=F 
37 
LOGICAL :: Flag_Opt_Geom 
38  
39 
! ====================================================================== 
40  
41 
LOGICAL :: debug 
42 
REAL(KREAL), ALLOCATABLE :: GradTmp(:) 
43 
REAL(KREAL), ALLOCATABLE :: GradCart(:) 
44 
REAL(KREAL), ALLOCATABLE :: x(:), y(:), z(:) 
45 
REAL(KREAL) :: Pi 
46  
47 
INTEGER(KINT) :: iat, i, j, IBeg,Idx 
48  
49 
REAL(KREAL), ALLOCATABLE :: XPrimRef(:), XPrim(:) ! NPrim 
50  
51 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
52 
! 
53 
!!!!!!!! PFL Debug 
54 
! 
55 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
56 
REAL(KREAL), ALLOCATABLE :: GeomTmp(:,:) !3,Nat 
57 
LOGICAL :: FALOBBT,FALOBPrimt,FAloBBTInv 
58 
LOGICAL :: DebugPFL 
59  
60 
! ====================================================================== 
61 
INTERFACE 
62 
function valid(string) result (isValid) 
63 
CHARACTER(*), intent(in) :: string 
64 
logical :: isValid 
65 
END function VALID 
66 
END INTERFACE 
67 
! ====================================================================== 
68  
69 
Pi=dacos(1.0d0) 
70  
71 
if (present(FOptGeom)) THEN 
72 
Flag_Opt_Geom=FOptGeom 
73 
ELSE 
74 
Flag_Opt_Geom=.FALSE. 
75 
END IF 
76  
77 
debug=valid('EGRAD') 
78 
debugPFL=valid('bakerPFL') 
79 
if (debug) Call Header("Entering Egrad") 
80  
81 
if (debug) Print *, 'EGrad.f90, L73, IOpt=', IOpt, ' Flag_Opt_Geom=',Flag_Opt_Geom 
82  
83 
Grad=0. 
84 
E=0. 
85  
86  
87 
ALLOCATE(GradCart(3*Nat)) 
88 
ALLOCATE(x(Nat),y(Nat),z(Nat)) 
89 
ALLOCATE(XPrimRef(NPrim),XPrim(NPrim)) 
90  
91 
SELECT CASE (COORD) 
92 
CASE ('ZMAT') 
93 
! In order to avoid problem with numbering and linear angles/diehedral, we convert the 
94 
! zmat into cartesian coordinates 
95 
! A remplacer par Int2Cart :) 
96 
Call Int2Cart(nat,indzmat,Geom,GeomCart) 
97 
CASE ('BAKER') 
98 
XPrimRef=XPrimitiveF(IGeom,:) 
99 
IF (Flag_Opt_Geom) THEN 
100 
IF (IOpt .EQ. 1) THEN 
101 
GeomCart(:,1) = XyzGeomI(IGeom,1,:) 
102 
GeomCart(:,2) = XyzGeomI(IGeom,2,:) 
103 
GeomCart(:,3) = XyzGeomI(IGeom,3,:) 
104 
ELSE 
105 
if (debug) WRITE(*,*) 'Before Call ConvertBakerInternal_cart; L125, Egrad.f90' 
106 
if (.NOT.present(GeomOld)) THEN 
107 
WRITE(*,*) "ERROR: GeomOld not passed as an argument" 
108 
STOP 
109 
END IF 
110 
WRITE(*,*) 'GeomOld=' 
111 
WRITE(*,'(12(1X,F6.3))') GeomOld 
112 
WRITE(*,*) 'Geom=' 
113 
WRITE(*,'(12(1X,F6.3))') Geom 
114 
if (.NOT.present(GeomCart_Old)) THEN 
115 
WRITE(*,*) "ERROR: GeomCart_Old not passed as an argument" 
116 
STOP 
117 
END IF 
118  
119 
WRITE(*,*) 'GeomCart_old=' 
120 
WRITE(*,'(20(1X,F6.3))') GeomCart_old(1:Nat,1:3) 
121 
Call ConvertBakerInternal_cart(GeomOld,Geom,GeomCart_old(:,1), & 
122 
GeomCart_old(:,2),GeomCart_old(:,3),x,y,z,XPrim,XPrimRef) 
123 
GeomCart(:,1) = x(:) 
124 
GeomCart(:,2) = y(:) 
125 
GeomCart(:,3) = z(:) 
126 
END IF 
127 
ELSE 
128 
IF (IOpt .EQ. 0) THEN 
129 
GeomCart(:,1) = XyzGeomF(IGeom,1,:) 
130 
GeomCart(:,2) = XyzGeomF(IGeom,2,:) 
131 
GeomCart(:,3) = XyzGeomF(IGeom,3,:) 
132 
ELSE 
133 
! IntCoordF(NGeomF,NCoord), XyzGeomF(NGeomF,3,Nat) 
134 
! Geom has to be converted into x,y,z 
135 
! GeomOld_all and XyzGeomF are old geometry in internal and cartesian coords resp. 
136 
DO I=1,NCoord ! these variables are used in ConvertBakerInternal_cart() 
137 
BTransInv_local(I,:) = BTransInvF(IGeom,I,:) 
138 
UMat_local(:,I) = UMatF(IGeom,:,I) 
139 
END DO 
140 
if (debug) Print *, 'I am before Call ConvertBakerInternal_cart; L160, Egrad.f90' 
141 
WRITE(*,*) 'GeomOld=' 
142 
WRITE(*,'(12(1X,F6.3))') GeomOld_all(IGeom,:) 
143 
WRITE(*,*) 'Geom=' 
144 
WRITE(*,'(12(1X,F6.3))') Geom 
145 
WRITE(*,*) 'DBG Egrad L165 GeomCart_old=' 
146 
DO I=1,Nat 
147 
WRITE(*,'(3(1X,F12.8))') XyzGeomF(IGeom,:,I) 
148 
END DO 
149  
150 
Call ConvertBakerInternal_cart(GeomOld_all(IGeom,:),Geom,XyzGeomF(IGeom,1,:), & 
151 
XyzGeomF(IGeom,2,:),XyzGeomF(IGeom,3,:),x,y,z,XPrim,XPrimRef) 
152 
if (debug) Print *, 'I am after Call ConvertBakerInternal_cart; L111, Egrad.f90' 
153 
DO I=1,NCoord ! these variables are used in ConvertBakerInternal_cart() 
154 
BTransInvF(IGeom,I,:) = BTransInv_local(I,:) 
155 
END DO 
156 
! GeomCart is actually XyzCoord and has INTENT(OUT) in this subroutine. 
157 
GeomCart(:,1) = x(:) 
158 
GeomCart(:,2) = y(:) 
159 
GeomCart(:,3) = z(:) 
160 
END IF ! matches IF (IOpt .EQ. 0) THEN 
161 
END IF 
162 
CASE ('MIXED') 
163 
! write(*,*) "Coucou 4mixed" 
164 
Call Mixed2Cart(Nat,indzmat,geom,GeomCart) 
165 
CASE ('CART') 
166 
! write(*,*) "Coucou 4cart" 
167 
if (debug) WRITE(*,*) "Coord=CART... in egrad" 
168 
GeomCart=reshape(Geom,(/Nat,3/)) 
169 
DO Iat=1,Nat 
170 
WRITE(IOTMP,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)),Geom(3*Iat2:3*Iat) 
171 
END DO 
172 
CASE ('HYBRID') 
173 
! write(*,*) "Coucou 4hybrid" 
174 
GeomCart=reshape(Geom,(/Nat,3/)) 
175 
CASE DEFAULT 
176 
WRITE (*,*) "Coord=",TRIM(COORD)," not recognized. Stop" 
177 
STOP 
178 
END SELECT 
179 
!WRITE(*,*) Nat 
180 
!WRITE(*,*) 'GeomCart in Egrad_baker' 
181 
!DO I=1,Nat 
182 
! WRITE(*,'(1X,A,3(1X,F12.8))') AtName(I),GeomCart(I,1:3) 
183 
!END DO 
184  
185 
SELECT CASE (Prog) 
186 
CASE ('GAUSSIAN') 
187 
! we call the Gaussian routine. 
188 
! it will return the energy and the gradient in cartesian coordinates. 
189 
! Gradient calculated at GeomCart geometry and stored in GradCart, which has INTENT(OUT) 
190 
Call egrad_gaussian(E,GeomCart,GradCart) 
191 
CASE ('MOPAC') 
192 
! we call the Mopac routine. 
193 
! it will return the energy and the gradient in cartesian coordinates. 
194 
! Gradient calculated at GeomCart geometry and stored in GradCart, which has INTENT(OUT) 
195 
Call egrad_mopac(E,GeomCart,GradCart) 
196 
CASE ('VASP') 
197 
Call egrad_vasp(E,Geomcart,GradCart) 
198 
CASE ('TURBOMOLE') 
199 
Call egrad_turbomole(E,Geomcart,GradCart) 
200 
CASE ('EXT') 
201 
Call egrad_ext(E,Geomcart,GradCart) 
202 
CASE ('TEST') 
203 
Call egrad_test(Nat,E,Geomcart,GradCart) 
204 
CASE ('CHAMFRE') 
205 
Call egrad_chamfre(Nat,E,Geomcart,GradCart) 
206 
CASE ('LEPS') 
207 
Call egrad_LEPS(Nat,E,Geomcart,GradCart) 
208 
END SELECT 
209 
if (debug) THEN 
210 
WRITE(*,*) 'DBG EGRAD, GradCart read' 
211 
DO I=1,Nat 
212 
Iat=I 
213 
if (renum) Iat=orderInv(I) 
214 
WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)),GradCart(3*I2:3*I) 
215 
END DO 
216 
END IF 
217  
218 
! If (PROG=="VASP") STOP 
219  
220  
221 
! We have the gradient in Cartesian coordinates... which is ok for COORD=Cart, Hybrid 
222 
! but we have to convert it into Zmat if COORD=Zmat 
223 
SELECT CASE (COORD) 
224 
CASE ("ZMAT") 
225 
ALLOCATE(GradTmp(3*Nat)) 
226 
if (debug) WRITE(*,*) "DBG EGRAD Calling CalcBMat_int" 
227 
CALL CalcBMat_int (nat, reshape(GeomCart,(/3,Nat/),ORDER=(/2,1/)), indzmat, dzdc, massat,atome) 
228  
229 
if (debug) WRITE(*,*) "DBG EGRAD Calling ConvGrad_Cart2Int" 
230 
CALL ConvGrad_Cart2Int(Nat,dzdc,IndZmat,GradCart,GradTmp) 
231  
232 
if (debug) WRITE(*,*) "DBG EGRAD Storing Grad" 
233 
Grad(1)=GradTmp(4) 
234 
Grad(2)=GradTmp(7) 
235 
! We might have troubles whith angles that are not in the [0:pi] range because, 
236 
! when converted to catesian coord, they do appear in the [0:Pi] range to gaussian... 
237 
! so that the derivative is not good, and a multiplicative factor should be added... 
238 
! 
239 
! This in fact should be taken care of in B mat calculation... 
240 
! 
241 
IF ((Geom(3).LE.0).OR.(Geom(3).GE.Pi)) THEN 
242 
Grad(3)=1.0d0*GradTmp(8) 
243 
ELSE 
244 
Grad(3)=GradTmp(8) 
245 
END IF 
246 
Idx=4 
247 
DO I=4,Nat 
248 
Grad(Idx)=GradTmp(3*i2) 
249 
Grad(Idx+1)=GradTmp(3*i1) 
250 
Grad(Idx+2)=GradTmp(3*i) 
251 
Idx=Idx+3 
252 
END DO 
253 
DEALLOCATE(GradTmp) 
254 
CASE ('BAKER') 
255 
! We have the gradient in Cartesian coordinates... which is ok for COORD=Cart, Hybrid 
256 
! but we have to convert it into internal coordinates if COORD=Baker 
257 
!!!!!!!!!!!!!!!!!!!! 
258 
! 
259 
! PFL Debug 
260 
! 
261 
!!!!!!!!!!!!!!!!!!!! 
262 
if (debugPFL) THEN 
263 
WRITE(*,*) "DBG PFL Egrad_Baker... computing dzdc" 
264 
if (.not.ALLOCATED(indzmat)) THEN 
265 
ALLOCATE(indzmat(Nat,5)) 
266 
IndZmat=0 
267 
Indzmat(1,1)=1 
268 
IndZmat(2,1)=2 
269 
IndZmat(2,2)=1 
270 
IndZmat(3,1)=3 
271 
IndZmat(3,2)=2 
272 
IndZmat(3,3)=1 
273 
IndZmat(4,1)=4 
274 
IndZmat(4,2)=3 
275 
IndZmat(4,3)=2 
276 
IndZmat(4,4)=1 
277 
END IF 
278 
IF (.NOT.ALLOCATED(DzDc)) THEN 
279 
ALLOCATE(DzDc(3,nat,3,Nat)) 
280 
END IF 
281 
if (.NOT.Allocated(GeomTmp)) Allocate(GeomTmp(3,Nat)) 
282 
DO I=1,Nat 
283 
GeomTmp(1,I)=GeomCart(OrderInv(I),1) 
284 
GeomTmp(2,I)=GeomCart(OrderInv(i),2) 
285 
GeomTmp(3,I)=GeomCart(OrderInv(i),3) 
286 
END DO 
287  
288 
DzDc=0.d0 
289 
CALL CalcBMat_int (nat, GeomTmp, indzmat, dzdc, massat,atome) 
290  
291 
END IF ! if (debugPFL) 
292 
!!!!!!!!!!!!!!!!!!!!!!!!! 
293 
! Debugging PFL 
294 
!!!!!!!!!!!!!!!!!!!!!!!! 
295  
296 
WRITE(*,*) "BTransInv_local Trans that is used originally" 
297 
DO J=1,3*Nat 
298 
WRITE(*,'(50(1X,F12.8))') BTransInv_local(:,J) 
299 
END DO 
300  
301 
WRITE(*,*) "Calculating actual values using GeomCart" 
302 
if (.NOT.Allocated(GeomTmp)) Allocate(GeomTmp(3,Nat)) 
303 
GeomTmp(1,:)=GeomCart(1:Nat,1) 
304 
GeomTmp(2,:)=GeomCart(1:Nat,2) 
305 
GeomTmp(3,:)=GeomCart(1:Nat,3) 
306  
307 
! Computing B^prim: 
308 
FAloBBT=.NOT.ALLOCATED(BBT) 
309 
FAloBBTInv=.NOT.ALLOCATED(BBT_inv) 
310 
FAloBPrimT=.NOT.ALLOCATED(BprimT) 
311 
if (FAloBPrimT) ALLOCATE(Bprimt(3*Nat,NPrim)) 
312 
if (FAloBBT) ALLOCATE(BBt(NCoord,NCoord)) 
313 
if (FAloBBTInv) ALLOCATE(BBt_inv(NCoord,NCoord)) 
314 
BprimT=0.d0 
315 
ScanCoord=>Coordinate 
316 
I=0 
317 
DO WHILE (Associated(ScanCoord%next)) 
318 
I=I+1 
319 
SELECT CASE (ScanCoord%Type) 
320 
CASE ('BOND') 
321 
CALL CONSTRAINTS_BONDLENGTH_DER(Nat,ScanCoord%at1,ScanCoord%AT2,GeomTmp,BprimT(1,I)) 
322 
CASE ('ANGLE') 
323 
CALL CONSTRAINTS_BONDANGLE_DER(Nat,ScanCoord%At1,ScanCoord%AT2,ScanCoord%At3,GeomTmp,BprimT(1,I)) 
324 
CASE ('DIHEDRAL') 
325 
CALL CONSTRAINTS_TORSION_DER2(Nat,ScanCoord%At1,ScanCoord%AT2,ScanCoord%At3,ScanCoord%At4,GeomTmp,BprimT(1,I)) 
326 
END SELECT 
327 
ScanCoord => ScanCoord%next 
328 
END DO 
329  
330 
if (debug) THEN 
331 
WRITE(*,*) "Bprim " 
332 
DO J=1,3*Nat 
333 
WRITE(*,'(50(1X,F12.6))') BprimT(J,:) 
334 
END DO 
335 
END IF 
336  
337 
BMat_BakerT = 0.d0 
338 
DO I=1,NCoord 
339 
DO J=1,NPrim 
340 
!BprimT is transpose of B^prim. 
341 
!B = UMat^T B^prim, B^T = (B^prim)^T UMat 
342 
BMat_BakerT(:,I)=BMat_BakerT(:,I)+BprimT(:,J)*UMat_local(J,I) 
343 
END DO 
344 
END DO 
345  
346 
BBT = 0.d0 
347 
DO I=1,NCoord 
348 
DO J=1,3*Nat ! BBT(:,I) forms BB^T 
349 
BBT(:,I) = BBT(:,I) + BMat_BakerT(J,:)*BMat_BakerT(J,I) 
350 
END DO 
351 
END DO 
352  
353 
BBT_inv = 0.d0 
354  
355 
Call GenInv(NCoord,BBT,BBT_inv,NCoord) 
356  
357 
! ALLOCATE(V(NCoord,NCoord)) 
358 
! tol = 1e10 
359 
! ! BBT is destroyed by GINVSE 
360 
! Call GINVSE(BBT,NCoord,NCoord,BBT_inv,NCoord,NCoord,NCoord,tol,V,NCoord,FAIL,IOOUT) 
361 
! DEALLOCATE(V) 
362 
! IF(Fail) Then 
363 
! WRITE (*,*) "Failed in generalized inverse, L220, ConvertBaker...f90" 
364 
! STOP 
365 
! END IF 
366  
367 
!Print *, 'BBT_inv=' 
368 
DO J=1,NCoord 
369 
! WRITE(*,'(50(1X,F10.2))') BBT_inv(:,J) 
370 
! Print *, BBT_inv(:,J) 
371 
END DO 
372 
!Print *, 'End of BBT_inv' 
373  
374 
! Calculation of (B^T)^1 = (BB^T)^1B: 
375 
BTransInv_local = 0.d0 
376 
DO I=1,3*Nat 
377 
DO J=1,NCoord 
378 
BTransInv_local(:,I)=BTransInv_local(:,I)+BBT_inv(:,J)*BMat_BakerT(I,J) 
379 
END DO 
380 
END DO 
381  
382 
Print *, 'BMat_BakerT=' 
383 
DO J=1,3*Nat 
384 
WRITE(*,'(50(1X,F12.8))') BMat_BakerT(:,J) 
385 
END DO 
386  
387 
WRITE(*,*) "BTransInv_local Trans corresponding to GeomCart" 
388 
DO J=1,3*Nat 
389 
WRITE(*,'(50(1X,F12.8))') BTransInv_local(:,J) 
390 
END DO 
391  
392 
if (FAloBPrimT) DEALLOCATE(Bprimt) 
393 
if (FAloBBT) DEALLOCATE(BBt) 
394 
if (FAloBBTInv) DEALLOCATE(BBt_inv) 
395  
396 
Grad=0.d0 
397 
DO I=1, 3*Nat 
398 
! Here Grad is gradient in Baker coordinates. GradCart(3*Nat) 
399 
Grad(:) = Grad(:) + BTransInv_local(:,I)*GradCart(I) !BTransInv_local = BTransInv in Opt_Geom.f90 
400 
END DO 
401 
!WRITE(*,*) 'In Egrad.f90, L263, Gradient:' 
402 
!WRITE(*,'(12(1X,F8.4))') Grad(1:NCoord) 
403  
404  
405 
CASE ('MIXED') 
406 
ALLOCATE(GradTmp(3*Nat)) 
407 
if (debug) WRITE(*,*) "DBG EGRAD Calling CalcBMat_mixed" 
408 
CALL CalcBMat_mixed (nat, reshape(GeomCart,(/3,Nat/),ORDER=(/2,1/)), indzmat, dzdc, massat,atome) 
409  
410 
if (Debug) THEN 
411 
WRITE(*,*) "DzDc" 
412 
DO I=1,Nat 
413 
DO J=1,3 
414 
WRITE(*,'(50(1X,F12.6))') dzdc(J,I,:,:) 
415 
END DO 
416 
END DO 
417 
END IF 
418 

419  
420 
CALL ConvGrad_Cart2Int(Nat,dzdc,IndZmat,GradCart,GradTmp) 
421 
DO I=1,Nat 
422 
WRITE(*,'(1X,3(1X,F10.5))') GradTmp(3*I2:3*I) 
423 
END DO 
424 
SELECT CASE (NCART) 
425 
CASE (1) 
426 
Grad(1:3)=GradTmp(1:3) 
427 
Grad(4)=GradTmp(4) 
428 
Grad(5)=GradTmp(7) 
429 
Grad(6)=GradTmp(8) 
430 
Idx=7 
431 
IBeg=4 
432 
CASE(2) 
433 
Grad(1:3)=GradTmp(1:3) 
434 
Grad(4:6)=GradTmp(4:6) 
435 
Grad(7)=GradTmp(7) 
436 
Grad(8)=GradTmp(8) 
437 
Idx=9 
438 
IBeg=4 
439 
CASE DEFAULT 
440 
Idx=1 
441 
IBeg=1 
442 
END SELECT 
443 
DO I=IBeg,Nat 
444 
Grad(Idx)=GradTmp(3*i2) 
445 
Grad(Idx+1)=GradTmp(3*i1) 
446 
Grad(Idx+2)=GradTmp(3*i) 
447 
Idx=Idx+3 
448 
END DO 
449 
DEALLOCATE(GradTmp) 
450 
CASE ("CART","HYBRID") 
451 
Grad=GradCart 
452 
CASE DEFAULT 
453 
WRITE(*,*) "Coord=",AdjustL(Trim(COORD))," not yet implemented in Egradients. STOP" 
454 
STOP 
455 
END SELECT 
456  
457 
DEALLOCATE(GradCart) 
458 
DEALLOCATE(x,y,z) 
459  
460 
IF (Debug) WRITE(*,*) 'DBG EGRAD grad',grad(1:NCoord) 
461 
if (debug) Call Header("Egrad Over") 
462  
463 
!999 CONTINUE 
464 
!if (.NOT.Ftmp) WRITE(*,*) 'We should not be here !!!!' 
465 
!STOP 
466 
! ====================================================================== 
467 
end subroutine Egrad 
468 