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

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