root / src / Egrad.f90 @ 9
History  View  Annotate  Download (15.3 kB)
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,Atome,massat, & 
8 
prog,NCart,XyzGeomF, XyzGeomI, renum, & 
9 
GeomOld_all,BTransInvF,BTransInv_local,UMatF,UMat_local & 
10 
, BprimT,ScanCoord, Coordinate,NPrim,BMat_BakerT,BBT,BBT_inv & 
11 
, 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  
167 
GeomCart=reshape(Geom,(/Nat,3/)) 
168 
if (debug) THEN 
169 
WRITE(*,*) "Coord=CART... in egrad" 
170 
DO Iat=1,Nat 
171 
WRITE(IOTMP,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)),Geom(3*Iat2:3*Iat) 
172 
END DO 
173 
END IF 
174 
CASE ('HYBRID') 
175 
! write(*,*) "Coucou 4hybrid" 
176 
GeomCart=reshape(Geom,(/Nat,3/)) 
177 
CASE DEFAULT 
178 
WRITE (*,*) "Coord=",TRIM(COORD)," not recognized. Stop" 
179 
STOP 
180 
END SELECT 
181 
!WRITE(*,*) Nat 
182 
!WRITE(*,*) 'GeomCart in Egrad_baker' 
183 
!DO I=1,Nat 
184 
! WRITE(*,'(1X,A,3(1X,F12.8))') AtName(I),GeomCart(I,1:3) 
185 
!END DO 
186  
187 
SELECT CASE (Prog) 
188 
CASE ('GAUSSIAN') 
189 
! we call the Gaussian routine. 
190 
! it will return the energy and the gradient in cartesian coordinates. 
191 
! Gradient calculated at GeomCart geometry and stored in GradCart, which has INTENT(OUT) 
192 
Call egrad_gaussian(E,GeomCart,GradCart) 
193 
CASE ('MOPAC') 
194 
! we call the Mopac routine. 
195 
! it will return the energy and the gradient in cartesian coordinates. 
196 
! Gradient calculated at GeomCart geometry and stored in GradCart, which has INTENT(OUT) 
197 
Call egrad_mopac(E,GeomCart,GradCart) 
198 
CASE ('VASP') 
199 
Call egrad_vasp(E,Geomcart,GradCart) 
200 
CASE ('SIESTA') 
201 
Call egrad_siesta(E,Geomcart,GradCart) 
202 
CASE ('TURBOMOLE') 
203 
Call egrad_turbomole(E,Geomcart,GradCart) 
204 
CASE ('EXT') 
205 
Call egrad_ext(E,Geomcart,GradCart) 
206 
CASE ('TEST') 
207 
Call egrad_test(Nat,E,Geomcart,GradCart) 
208 
CASE ('TEST2D') 
209 
Call egrad_test_2D(Nat,E,Geomcart,GradCart) 
210 
CASE ('CHAMFRE') 
211 
Call egrad_chamfre(Nat,E,Geomcart,GradCart) 
212 
CASE ('LEPS') 
213 
Call egrad_LEPS(Nat,E,Geomcart,GradCart) 
214 
END SELECT 
215 
if (debug) THEN 
216 
WRITE(*,*) 'DBG EGRAD, GradCart read' 
217 
DO I=1,Nat 
218 
Iat=I 
219 
if (renum) Iat=orderInv(I) 
220 
WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)),GradCart(3*I2:3*I) 
221 
END DO 
222 
END IF 
223  
224 
! If (PROG=="VASP") STOP 
225  
226  
227 
! We have the gradient in Cartesian coordinates... which is ok for COORD=Cart, Hybrid 
228 
! but we have to convert it into Zmat if COORD=Zmat 
229 
SELECT CASE (COORD) 
230 
CASE ("ZMAT") 
231 
ALLOCATE(GradTmp(3*Nat)) 
232 
if (debug) WRITE(*,*) "DBG EGRAD Calling CalcBMat_int" 
233 
CALL CalcBMat_int (nat, reshape(GeomCart,(/3,Nat/),ORDER=(/2,1/)), indzmat, dzdc, massat,atome) 
234  
235 
if (debug) WRITE(*,*) "DBG EGRAD Calling ConvGrad_Cart2Int" 
236 
CALL ConvGrad_Cart2Int(Nat,dzdc,IndZmat,GradCart,GradTmp) 
237  
238 
if (debug) WRITE(*,*) "DBG EGRAD Storing Grad" 
239 
Grad(1)=GradTmp(4) 
240 
Grad(2)=GradTmp(7) 
241 
! We might have troubles whith angles that are not in the [0:pi] range because, 
242 
! when converted to catesian coord, they do appear in the [0:Pi] range to gaussian... 
243 
! so that the derivative is not good, and a multiplicative factor should be added... 
244 
! 
245 
! This in fact should be taken care of in B mat calculation... 
246 
! 
247 
IF ((Geom(3).LE.0).OR.(Geom(3).GE.Pi)) THEN 
248 
Grad(3)=1.0d0*GradTmp(8) 
249 
ELSE 
250 
Grad(3)=GradTmp(8) 
251 
END IF 
252 
Idx=4 
253 
DO I=4,Nat 
254 
Grad(Idx)=GradTmp(3*i2) 
255 
Grad(Idx+1)=GradTmp(3*i1) 
256 
Grad(Idx+2)=GradTmp(3*i) 
257 
Idx=Idx+3 
258 
END DO 
259 
DEALLOCATE(GradTmp) 
260 
CASE ('BAKER') 
261 
! We have the gradient in Cartesian coordinates... which is ok for COORD=Cart, Hybrid 
262 
! but we have to convert it into internal coordinates if COORD=Baker 
263 
!!!!!!!!!!!!!!!!!!!! 
264 
! 
265 
! PFL Debug 
266 
! 
267 
!!!!!!!!!!!!!!!!!!!! 
268 
if (debugPFL) THEN 
269 
WRITE(*,*) "DBG PFL Egrad_Baker... computing dzdc" 
270 
if (.not.ALLOCATED(indzmat)) THEN 
271 
ALLOCATE(indzmat(Nat,5)) 
272 
IndZmat=0 
273 
Indzmat(1,1)=1 
274 
IndZmat(2,1)=2 
275 
IndZmat(2,2)=1 
276 
IndZmat(3,1)=3 
277 
IndZmat(3,2)=2 
278 
IndZmat(3,3)=1 
279 
IndZmat(4,1)=4 
280 
IndZmat(4,2)=3 
281 
IndZmat(4,3)=2 
282 
IndZmat(4,4)=1 
283 
END IF 
284 
IF (.NOT.ALLOCATED(DzDc)) THEN 
285 
ALLOCATE(DzDc(3,nat,3,Nat)) 
286 
END IF 
287 
if (.NOT.Allocated(GeomTmp)) Allocate(GeomTmp(3,Nat)) 
288 
DO I=1,Nat 
289 
GeomTmp(1,I)=GeomCart(OrderInv(I),1) 
290 
GeomTmp(2,I)=GeomCart(OrderInv(i),2) 
291 
GeomTmp(3,I)=GeomCart(OrderInv(i),3) 
292 
END DO 
293  
294 
DzDc=0.d0 
295 
CALL CalcBMat_int (nat, GeomTmp, indzmat, dzdc, massat,atome) 
296  
297 
END IF ! if (debugPFL) 
298 
!!!!!!!!!!!!!!!!!!!!!!!!! 
299 
! Debugging PFL 
300 
!!!!!!!!!!!!!!!!!!!!!!!! 
301  
302 
WRITE(*,*) "BTransInv_local Trans that is used originally" 
303 
DO J=1,3*Nat 
304 
WRITE(*,'(50(1X,F12.8))') BTransInv_local(:,J) 
305 
END DO 
306  
307 
WRITE(*,*) "Calculating actual values using GeomCart" 
308 
if (.NOT.Allocated(GeomTmp)) Allocate(GeomTmp(3,Nat)) 
309 
GeomTmp(1,:)=GeomCart(1:Nat,1) 
310 
GeomTmp(2,:)=GeomCart(1:Nat,2) 
311 
GeomTmp(3,:)=GeomCart(1:Nat,3) 
312  
313 
! Computing B^prim: 
314 
FAloBBT=.NOT.ALLOCATED(BBT) 
315 
FAloBBTInv=.NOT.ALLOCATED(BBT_inv) 
316 
FAloBPrimT=.NOT.ALLOCATED(BprimT) 
317 
if (FAloBPrimT) ALLOCATE(Bprimt(3*Nat,NPrim)) 
318 
if (FAloBBT) ALLOCATE(BBt(NCoord,NCoord)) 
319 
if (FAloBBTInv) ALLOCATE(BBt_inv(NCoord,NCoord)) 
320 
BprimT=0.d0 
321 
ScanCoord=>Coordinate 
322 
I=0 
323 
DO WHILE (Associated(ScanCoord%next)) 
324 
I=I+1 
325 
SELECT CASE (ScanCoord%Type) 
326 
CASE ('BOND') 
327 
CALL CONSTRAINTS_BONDLENGTH_DER(Nat,ScanCoord%at1,ScanCoord%AT2,GeomTmp,BprimT(1,I)) 
328 
CASE ('ANGLE') 
329 
CALL CONSTRAINTS_BONDANGLE_DER(Nat,ScanCoord%At1,ScanCoord%AT2,ScanCoord%At3,GeomTmp,BprimT(1,I)) 
330 
CASE ('DIHEDRAL') 
331 
CALL CONSTRAINTS_TORSION_DER2(Nat,ScanCoord%At1,ScanCoord%AT2,ScanCoord%At3,ScanCoord%At4,GeomTmp,BprimT(1,I)) 
332 
END SELECT 
333 
ScanCoord => ScanCoord%next 
334 
END DO 
335  
336 
if (debug) THEN 
337 
WRITE(*,*) "Bprim " 
338 
DO J=1,3*Nat 
339 
WRITE(*,'(50(1X,F12.6))') BprimT(J,:) 
340 
END DO 
341 
END IF 
342  
343 
BMat_BakerT = 0.d0 
344 
DO I=1,NCoord 
345 
DO J=1,NPrim 
346 
!BprimT is transpose of B^prim. 
347 
!B = UMat^T B^prim, B^T = (B^prim)^T UMat 
348 
BMat_BakerT(:,I)=BMat_BakerT(:,I)+BprimT(:,J)*UMat_local(J,I) 
349 
END DO 
350 
END DO 
351  
352 
BBT = 0.d0 
353 
DO I=1,NCoord 
354 
DO J=1,3*Nat ! BBT(:,I) forms BB^T 
355 
BBT(:,I) = BBT(:,I) + BMat_BakerT(J,:)*BMat_BakerT(J,I) 
356 
END DO 
357 
END DO 
358  
359 
BBT_inv = 0.d0 
360  
361 
Call GenInv(NCoord,BBT,BBT_inv,NCoord) 
362  
363 
! ALLOCATE(V(NCoord,NCoord)) 
364 
! tol = 1e10 
365 
! ! BBT is destroyed by GINVSE 
366 
! Call GINVSE(BBT,NCoord,NCoord,BBT_inv,NCoord,NCoord,NCoord,tol,V,NCoord,FAIL,IOOUT) 
367 
! DEALLOCATE(V) 
368 
! IF(Fail) Then 
369 
! WRITE (*,*) "Failed in generalized inverse, L220, ConvertBaker...f90" 
370 
! STOP 
371 
! END IF 
372  
373 
!Print *, 'BBT_inv=' 
374 
DO J=1,NCoord 
375 
! WRITE(*,'(50(1X,F10.2))') BBT_inv(:,J) 
376 
! Print *, BBT_inv(:,J) 
377 
END DO 
378 
!Print *, 'End of BBT_inv' 
379  
380 
! Calculation of (B^T)^1 = (BB^T)^1B: 
381 
BTransInv_local = 0.d0 
382 
DO I=1,3*Nat 
383 
DO J=1,NCoord 
384 
BTransInv_local(:,I)=BTransInv_local(:,I)+BBT_inv(:,J)*BMat_BakerT(I,J) 
385 
END DO 
386 
END DO 
387  
388 
Print *, 'BMat_BakerT=' 
389 
DO J=1,3*Nat 
390 
WRITE(*,'(50(1X,F12.8))') BMat_BakerT(:,J) 
391 
END DO 
392  
393 
WRITE(*,*) "BTransInv_local Trans corresponding to GeomCart" 
394 
DO J=1,3*Nat 
395 
WRITE(*,'(50(1X,F12.8))') BTransInv_local(:,J) 
396 
END DO 
397  
398 
if (FAloBPrimT) DEALLOCATE(Bprimt) 
399 
if (FAloBBT) DEALLOCATE(BBt) 
400 
if (FAloBBTInv) DEALLOCATE(BBt_inv) 
401  
402 
Grad=0.d0 
403 
DO I=1, 3*Nat 
404 
! Here Grad is gradient in Baker coordinates. GradCart(3*Nat) 
405 
Grad(:) = Grad(:) + BTransInv_local(:,I)*GradCart(I) !BTransInv_local = BTransInv in Opt_Geom.f90 
406 
END DO 
407 
!WRITE(*,*) 'In Egrad.f90, L263, Gradient:' 
408 
!WRITE(*,'(12(1X,F8.4))') Grad(1:NCoord) 
409  
410  
411 
CASE ('MIXED') 
412 
ALLOCATE(GradTmp(3*Nat)) 
413 
if (debug) WRITE(*,*) "DBG EGRAD Calling CalcBMat_mixed" 
414 
CALL CalcBMat_mixed (nat, reshape(GeomCart,(/3,Nat/),ORDER=(/2,1/)), indzmat, dzdc, massat,atome) 
415  
416 
if (Debug) THEN 
417 
WRITE(*,*) "DzDc" 
418 
DO I=1,Nat 
419 
DO J=1,3 
420 
WRITE(*,'(50(1X,F12.6))') dzdc(J,I,:,:) 
421 
END DO 
422 
END DO 
423 
END IF 
424 

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