1*
2* $Id$
3*
4
5*    ************************************
6*    *					*
7*    *	    gen_revPBE_unrestricted	*
8*    *					*
9*    ************************************
10*
11*    This function returns the revPBE exchange-correlation
12*  energy density, xce, and its derivatives with respect
13*  to nup, ndn, |grad nup|, |grad ndn|, and |grad n|.
14*
15*   Entry - n2ft3d     : number of grid points
16*           dn_in(*,2) : spin densites nup and ndn
17*           agr_in(*,3): |grad nup|, |grad ndn|, and |grad n|
18*           x_parameter: scale parameter for exchange
19*           c_parameter: scale parameter for correlation
20*
21*   Exit - xce(*)  : revPBE energy density
22*        - fn(*,2) : d(n*xce)/dnup, d(n*xce)/dndn
23*        - fdn(*,3): d(n*xce)/d|grad nup|, d(n*xce)/d|grad ndn|
24*                    d(n*xce)/d|grad n|
25
26      subroutine gen_revPBE_BW_unrestricted(n2ft3d,
27     >                           dn_in,agr_in,
28     >                           x_parameter,c_parameter,
29     >                           xce,fn,fdn)
30      implicit none
31
32      integer n2ft3d
33      real*8 dn_in(n2ft3d,2)
34      real*8 agr_in(n2ft3d,3)
35      real*8 x_parameter,c_parameter
36      real*8 xce(n2ft3d)
37      real*8 fn(n2ft3d,2)
38      real*8 fdn(n2ft3d,3)
39
40*     **** Density cutoff parameter ****
41      real*8 DNS_CUT,ETA,ETA2,alpha_zeta,alpha_zeta2
42      parameter (DNS_CUT = 1.0d-20)
43      parameter (ETA=1.0d-20)
44      parameter (ETA2=1.0d-14)
45      parameter (alpha_zeta=(1.0d0-ETA2))
46      parameter (alpha_zeta2=(1.0d0-ETA2))
47
48c     ***** revPBE GGA exchange constants ******
49      real*8 MU,KAPPA
50      parameter (MU    = 0.2195149727645171d0)
51      !parameter (KAPPA = 0.8040000000000000d0)
52      parameter (KAPPA = 1.2450000000000000d0)
53
54c     ****** PBE96 GGA correlation constants ******
55      real*8 GAMMA,BETA,BOG
56      parameter (GAMMA	= 0.031090690869655d0)
57      parameter (BETA	= 0.066724550603149d0)
58      !parameter (BETA	= 0.066725d0)
59      parameter (BOG    = BETA/GAMMA)
60
61
62c     ****** Perdew-Wang92 LDA correlation coefficients *******
63      real*8 GAM,iGAM,FZZ,iFZZ
64      parameter (GAM  	= 0.519842099789746329d0)
65      parameter (iGAM  	= 1.0d0/GAM)
66      parameter (FZZ    = (8.0d0/(9.0d0*GAM)) )
67      parameter (iFZZ    = 0.125d0*9.0d0*GAM)
68
69      real*8 A_1,A1_1,B1_1,B2_1,B3_1,B4_1
70      parameter (A_1  = 0.0310907d0)
71      !parameter (A_1  = 0.031091d0)
72      parameter (A1_1 =	0.2137000d0)
73      parameter (B1_1 =	7.5957000d0)
74      parameter (B2_1 =	3.5876000d0)
75      parameter (B3_1 =	1.6382000d0)
76      parameter (B4_1 =	0.4929400d0)
77
78      real*8 A_2,A1_2,B1_2,B2_2,B3_2,B4_2
79      parameter (A_2  =  0.01554535d0)
80      !parameter (A_2  =  0.015545d0)
81      parameter (A1_2 =	 0.20548000d0)
82      parameter (B1_2 =	14.11890000d0)
83      parameter (B2_2 =	 6.19770000d0)
84      parameter (B3_2 =	 3.36620000d0)
85      parameter (B4_2 =	 0.62517000d0)
86
87      real*8 A_3,A1_3,B1_3,B2_3,B3_3,B4_3
88      parameter (A_3  =  0.0168869d0)
89      !parameter (A_3  =  0.016887d0)
90      parameter (A1_3 =	 0.1112500d0)
91      parameter (B1_3 =	10.3570000d0)
92      parameter (B2_3 =	 3.6231000d0)
93      parameter (B3_3 =	 0.8802600d0)
94      parameter (B4_3 =	 0.4967100d0)
95
96c     **** other constants ****
97      real*8 onethird,fourthird,fivethird,onesixthm
98      real*8 twothird,sevensixthm
99      real*8 onethirdm
100      parameter (onethird=1.0d0/3.0d0)
101      parameter (onethirdm=-1.0d0/3.0d0)
102      parameter (twothird=2.0d0/3.0d0)
103      parameter (fourthird=4.0d0/3.0d0)
104      parameter (fivethird=5.0d0/3.0d0)
105      parameter (onesixthm=-1.0d0/6.0d0)
106      parameter (sevensixthm=-7.0d0/6.0d0)
107
108c     **** local variables ****
109      integer i
110      real*8 n,agr
111      real*8 nup,agrup
112      real*8 ndn,agrdn
113      real*8 kf,ks,s,P0,n_onethird,pi,rs_scale
114      real*8 rs         ! Wigner radius
115      real*8 rss        ! rss  = sqrt(rs)
116      real*8 rs_n       ! rs_n = n*drs/dn
117      real*8 t,t2,t4,t6
118      real*8 t_nup      ! t_nup = n*dt/dnup
119      real*8 t_ndn      ! t_ndn = n*dt/dndn
120      real*8 t_agr      ! t_agr = n*dt/dagr
121      real*8 zet,twoksg
122      real*8 zet_nup    ! zet_nup = n*dzet/dnup
123      real*8 zet_ndn    ! zet_nup = n*dzet/dnup
124      real*8 zetp_1_3,zetm_1_3
125      real*8 zetpm_1_3,zetmm_1_3
126      real*8 phi,phi3,phi4
127      real*8 phi_zet
128      real*8 A,A2
129      real*8 A_phi,A_ec_lda
130      real*8 Q0,Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8
131      real*8 PON,FZ,z4
132      real*8 tau
133      real*8 F
134      real*8 Fs                ! dF/ds
135      real*8 Hpbe
136      real*8 Hpbe_t            ! dHpbe/dt
137      real*8 Hpbe_phi          ! dHpbe/dphi
138      real*8 Hpbe_ec_lda       ! dHpbe/d(ec_lda)
139      real*8 Hpbe_nup,Hpbe_ndn ! n*dHpbe/dnup, n*dHpbe/dndn
140      real*8 Ipbe
141      real*8 Ipbe_t,Ipbe_A     ! dIpbe/dt, dIpbe/dA
142
143      real*8 exup,exdn,ex,ex_lda
144      real*8 ecu,ecp,eca,ec,ec_lda
145      real*8 ecu_rs,ecp_rs,eca_rs
146      real*8 ec_lda_rs,ec_lda_zet  ! d(ec_lda)/drs, d(ec_lda)/dzet
147      real*8 ec_lda_nup,ec_lda_ndn ! n*d(ec_lda)/dnup, n*d(ec_lda)/dndn
148      real*8 fnxup,fdnxup          ! d(n*ex)/dnup, d(n*ex)/dndn
149      real*8 fnxdn,fdnxdn          ! d(n*ex)/d|grad nup|, d(n*ex)/d|grad ndn|
150      real*8 fncup,fncdn           ! d(n*ec)/dnup, d(n*ec)/dndn
151      real*8 fdnx_const
152
153      pi = 4.0d0*datan(1.0d0)
154      rs_scale = (0.75d0/pi)**onethird
155      fdnx_const = -3.0d0/(8.0d0*pi)
156
157
158!$OMP DO
159      do i=1,n2ft3d
160         nup     = dn_in(i,1)+ETA
161         agrup   = agr_in(i,1)
162
163         ndn     = dn_in(i,2)+ETA
164         agrdn   = agr_in(i,2)
165
166c        ****************************************************************
167c        ***** calculate polarized Exchange energies and potentials *****
168c        ****************************************************************
169
170c        ************
171c        **** up ****
172c        ************
173         n     = 2.0d0*nup
174         agr   = 2.0d0*agrup
175
176         n_onethird = (3.0d0*n/pi)**onethird
177         ex_lda     = -0.75d0*n_onethird
178
179         kf = (3.0d0*pi*pi*n)**onethird
180         s  = agr/(2.0d0*kf*n)
181         P0 = 1.0d0 + (MU/KAPPA)*s*s
182
183         F   = (1.0d0 + KAPPA - KAPPA/P0)
184         Fs  = 2.0d0*MU/(P0*P0)*s
185
186         exup = ex_lda*F
187         fnxup = fourthird*(exup - ex_lda*Fs*s)
188         fdnxup = fdnx_const*Fs
189
190c        **************
191c        **** down ****
192c        **************
193         n     = 2.0d0*ndn
194         agr   = 2.0d0*agrdn
195
196         n_onethird = (3.0d0*n/pi)**onethird
197         ex_lda     = -0.75d0*n_onethird
198
199         kf = (3.0d0*pi*pi*n)**onethird
200         s  = agr/(2.0d0*kf*n)
201         P0 = 1.0d0 + (MU/KAPPA)*s*s
202
203         F   = (1.0d0 + KAPPA - KAPPA/P0)
204         Fs  = 2.0d0*MU/(P0*P0)*s
205
206         exdn   = ex_lda*F
207         fnxdn  = fourthird*(exdn - ex_lda*Fs*s)
208         fdnxdn = fdnx_const*Fs
209
210         n = nup+ndn
211
212         ex = (exup*nup+ exdn*ndn)/ n
213
214c        *******************************************************************
215c        ***** calculate polarized correlation energies and potentials *****
216c        *******************************************************************
217         agr   = agr_in(i,3)
218
219         zet = (nup-ndn)/n
220c         if (zet.gt.0.0d0) zet = zet - ETA2
221c         if (zet.lt.0.0d0) zet = zet + ETA2
222c        if (dabs(dn_in(i,2)).gt.DNS_CUT) zet_nup =  2*ndn/n**2
223c        if (dabs(dn_in(i,1)).gt.DNS_CUT) zet_ndn = -2*nup/n**2
224c        if (dabs(dn_in(i,2)).gt.DNS_CUT) zet_nup =  2*ndn/n
225c        zet_nup =  2*ndn/n
226c        zet_ndn = -2*nup/n
227         zet_nup = -(zet - 1.0d0)
228         zet_ndn = -(zet + 1.0d0)
229         zetpm_1_3 = (1.0d0+zet*alpha_zeta)**onethirdm
230         zetmm_1_3 = (1.0d0-zet*alpha_zeta)**onethirdm
231         zetp_1_3  = (1.0d0+zet*alpha_zeta)*zetpm_1_3**2
232         zetm_1_3  = (1.0d0-zet*alpha_zeta)*zetmm_1_3**2
233
234
235         phi = 0.5d0*( zetp_1_3**2 + zetm_1_3**2)
236         phi_zet = alpha_zeta*( zetpm_1_3 - zetmm_1_3)/3.0d0
237         F =(  (1.0d0+zet*alpha_zeta)*zetp_1_3
238     >       + (1.0d0-zet*alpha_zeta)*zetm_1_3
239     >       - 2.0d0)*iGAM
240
241         FZ = (zetp_1_3 - zetm_1_3)*(alpha_zeta*fourthird*iGAM)
242
243
244
245*        **** calculate Wigner radius ****
246         rs    = rs_scale/(n**onethird)
247         rss   = dsqrt(rs)
248
249*        **** calculate n*drs/dn ****
250c        rs_n = onethirdm*rs/n
251         rs_n = onethirdm*rs
252
253
254
255c        **** calculate t ****
256         kf = (3.0d0*pi*pi*n)**onethird
257         ks = dsqrt(4.0d0*kf/pi)
258
259         twoksg = 2.0d0*ks*phi
260
261         t  = agr/(twoksg*n)
262
263*        *** calculate n*dt/dnup, n*dt/dndn, n*dt/d|grad n| ****
264         t_nup = sevensixthm*t - (phi_zet)*(zet_nup)*t/phi
265         t_ndn = sevensixthm*t - (phi_zet)*(zet_ndn)*t/phi
266         t_agr  = 1.0d0/(twoksg)
267
268
269
270
271c        **************************************************
272c        ***** compute LSDA correlation energy density ****
273c        **************************************************
274         call LSDT2(A_1,A1_1,B1_1,B2_1,B3_1,B4_1,rss,ecu,ecu_rs)
275         call LSDT2(A_2,A1_2,B1_2,B2_2,B3_2,B4_2,rss,ecp,ecp_rs)
276         call LSDT2(A_3,A1_3,B1_3,B2_3,B3_3,B4_3,rss,eca,eca_rs)
277
278         z4 = zet**4
279
280         ec_lda = ecu*(1.0d0-F*z4)
281     >          + ecp*F*z4
282     >          - eca*F*(1.0d0-z4)/FZZ
283
284         ec_lda_rs = ecu_rs*(1.0d0-F*z4)
285     >             + ecp_rs*F*z4
286     >             - eca_rs*F*(1.0d0-z4)/FZZ
287
288         ec_lda_zet = (4.0d0*(zet**3)*F + FZ*z4)*(ecp-ecu+eca*iFZZ)
289     >              - FZ*eca*iFZZ
290
291
292
293
294c        ********************************************
295c        **** calculate PBE96 correlation energy ****
296c        ********************************************
297         phi3 = phi**3
298         phi4 = phi3*phi
299         PON  = -ec_lda/(phi3*GAMMA)
300         tau  = DEXP(PON)
301
302         A = BOG/(tau-1.0d0+ETA)
303         A2 = A*A
304         t2 = t*t
305         t4 = t2*t2
306         t6 = t4*t2
307         Q4 = 1.0d0 + A*t2
308         Q5 = 1.0d0 + 2.0d0*A*t2
309         Q6 = 2.0d0 + A*t2
310         Q7 = 1.0d0+A*t2+A2*t4
311         Q8 = Q7*Q7
312
313         Ipbe = 1.0d0 + BOG*t2*Q4/Q7
314         Hpbe = GAMMA*phi3*DLOG(Ipbe)
315
316         Ipbe_t =  BOG*(2.0d0*t)*Q5/Q8
317         Ipbe_A = -BOG*(A*t6)   *Q6/Q8
318
319         A_ec_lda  = tau/(BETA*phi3)*A2
320         A_phi     = -3.0d0*ec_lda*tau/(BETA*phi4)*A2
321
322
323         Hpbe_ec_lda = (GAMMA*phi3/Ipbe)*Ipbe_A*A_ec_lda
324
325         Hpbe_phi    = 3.0d0*Hpbe/phi
326     >               + (GAMMA*phi3/Ipbe)*Ipbe_A*A_phi
327
328         Hpbe_t      = (GAMMA*phi3/Ipbe)*Ipbe_t
329
330         ec_lda_nup = ec_lda_zet
331     >              - zet * ec_lda_zet
332     >              + rs_n * ec_lda_rs
333         ec_lda_ndn = -ec_lda_zet
334     >              - zet  * ec_lda_zet
335     >              + rs_n * ec_lda_rs
336
337
338
339         Hpbe_nup  = ec_lda_nup   * Hpbe_ec_lda
340     >          + phi_zet*zet_nup * Hpbe_phi
341     >          + t_nup           * Hpbe_t
342
343         Hpbe_ndn  = ec_lda_ndn   * Hpbe_ec_lda
344     >          + phi_zet*zet_ndn * Hpbe_phi
345     >          + t_ndn           * Hpbe_t
346
347
348
349         ec = ec_lda + Hpbe
350
351         fncup  = ec + (ec_lda_nup + Hpbe_nup)
352         fncdn  = ec + (ec_lda_ndn + Hpbe_ndn)
353
354         xce(i)   = x_parameter*ex     + c_parameter*ec
355         fn(i,1)  = x_parameter*fnxup  + c_parameter*fncup
356         fn(i,2)  = x_parameter*fnxdn  + c_parameter*fncdn
357
358         fdn(i,1) = x_parameter*fdnxup
359         fdn(i,2) = x_parameter*fdnxdn
360         fdn(i,3) = c_parameter*t_agr*Hpbe_t
361
362      end do
363!$OMP END DO
364
365
366
367      return
368      end
369
370
371
372
373
374*    ************************************
375*    *					*
376*    *	    gen_revPBE_BW_restricted	*
377*    *					*
378*    ************************************
379*
380*   This routine calculates the revPBE exchange-correlation
381*   potential(xcp) and energy density(xce).
382*
383*
384*   Entry - n2ft3d     : number of grid points
385*           rho_in(*) :  density (nup+ndn)
386*           agr_in(*): |grad rho_in|
387*           x_parameter: scale parameter for exchange
388*           c_parameter: scale parameter for correlation
389*
390*     Exit  - xce(n2ft3d) : revPBE exchange correlation energy density
391*             fn(n2ft3d)  : d(n*xce)/dn
392*             fdn(n2ft3d) : d(n*xce/d|grad n|
393*
394      subroutine gen_revPBE_BW_restricted(n2ft3d,rho_in,agr_in,
395     >                                x_parameter,c_parameter,
396     >                                xce,fn,fdn)
397      implicit none
398
399      integer    n2ft3d
400      real*8     rho_in(n2ft3d)
401      real*8     agr_in(n2ft3d)
402      real*8     x_parameter,c_parameter
403      real*8     xce(n2ft3d)
404      real*8     fn(n2ft3d)
405      real*8     fdn(n2ft3d)
406
407
408*     **** Density cutoff parameter ****
409      real*8 DNS_CUT,ETA
410      parameter (DNS_CUT = 1.0d-20)
411      parameter (ETA     = 1.0d-20)
412
413c     ***** revPBE GGA exchange constants ******
414      real*8 MU,KAPPA
415      parameter (MU    = 0.2195149727645171d0)
416      !parameter (KAPPA = 0.8040000000000000d0)
417      parameter (KAPPA = 1.2450000000000000d0)
418
419c     ****** revPBE GGA correlation constants ******
420      real*8 GAMMA,BETA,BOG
421      parameter (GAMMA	= 0.031090690869655d0)
422      parameter (BETA	= 0.066724550603149d0)
423      parameter (BOG    = BETA/GAMMA)
424
425
426c     ****** Perdew-Wang92 LDA correlation coefficients *******
427      real*8 A_1,A1_1,B1_1,B2_1,B3_1,B4_1
428      parameter (A_1  = 0.0310907d0)
429      parameter (A1_1 =	0.2137000d0)
430      parameter (B1_1 =	7.5957000d0)
431      parameter (B2_1 =	3.5876000d0)
432      parameter (B3_1 =	1.6382000d0)
433      parameter (B4_1 =	0.4929400d0)
434
435      real*8 A_2,A1_2,B1_2,B2_2,B3_2,B4_2
436      parameter (A_2  =  0.01554535d0)
437      parameter (A1_2 =	 0.20548000d0)
438      parameter (B1_2 =	14.11890000d0)
439      parameter (B2_2 =	 6.19770000d0)
440      parameter (B3_2 =	 3.36620000d0)
441      parameter (B4_2 =	 0.62517000d0)
442
443      real*8 A_3,A1_3,B1_3,B2_3,B3_3,B4_3
444      parameter (A_3  =  0.0168869d0)
445      parameter (A1_3 =	 0.1112500d0)
446      parameter (B1_3 =	10.3570000d0)
447      parameter (B2_3 =	 3.6231000d0)
448      parameter (B3_3 =	 0.8802600d0)
449      parameter (B4_3 =	 0.4967100d0)
450
451c     **** other constants ****
452      real*8 onethird,fourthird,sevensixths
453      parameter (onethird=1.0d0/3.0d0)
454      parameter (fourthird=4.0d0/3.0d0)
455      parameter (sevensixths=7.0d0/6.0d0)
456
457c     **** local variables ****
458      integer i
459      real*8 n,agr
460      real*8 kf,ks,s,P0,n_onethird,pi,rs_scale
461      real*8 fdnx_const
462      real*8 rs,rss,t,t2,t4,t6
463      real*8 Q0,Q1,Q2,Q3,Q4,Q5,Q8,Q9,B
464      real*8 Ht
465      real*8 B_ec,Hrs,H_B
466      real*8 F,Fs
467
468      real*8 ex_lda,ec_lda
469      real*8 ec_lda_rs
470      real*8 ex,ec,H
471      real*8 fnx,fdnx,fnc,fdnc
472
473
474      pi         = 4.0d0*datan(1.0d0)
475      rs_scale   = (0.75d0/pi)**onethird
476      fdnx_const = -3.0d0/(8.0d0*pi)
477
478!$OMP DO
479      do i=1,n2ft3d
480         n     = rho_in(i)+ETA
481         agr   = agr_in(i)
482
483c        ***** calculate unpolarized Exchange energies and potentials *****
484         n_onethird = (3.0d0*n/pi)**onethird
485         ex_lda     = -0.75d0*n_onethird
486
487         kf = (3.0d0*pi*pi*n)**onethird
488         s  = agr/(2.0d0*kf*n)
489         P0 = 1.0d0 + (MU/KAPPA)*s*s
490
491c        if (n.gt.DNS_CUT) then
492c           F   = (1.0d0 + KAPPA - KAPPA/P0)
493c           Fs  = 2.0d0*MU/(P0*P0)*s
494c        else
495c           F   = 1.0d0
496c           Fs  = 0.0d0
497c        end if
498         F   = (1.0d0 + KAPPA - KAPPA/P0)
499         Fs  = 2.0d0*MU/(P0*P0)*s
500
501         ex   = ex_lda*F
502         fnx  = fourthird*(ex - ex_lda*Fs*s)
503         fdnx = fdnx_const*Fs
504
505
506*        *********************************************************************
507c        ***** calculate unpolarized correlation energies and potentials *****
508*        *********************************************************************
509
510c        **** calculate rs and t ****
511         rs    = rs_scale/(n**onethird)
512         rss   = dsqrt(rs)
513
514         kf = (3.0d0*pi*pi*n)**onethird
515         ks = dsqrt(4.0d0*kf/pi)
516         t  = agr/(2.0d0*ks*n)
517
518
519c        **** unpolarized LDA correlation energy ****
520c        **** ec_p = correlation energy          ****
521c        ****   ec_p_rs = dec_p/drs              ****
522c        ****   uc_p    = dec_p/dn               ****
523         call LSDT2(A_1,A1_1,B1_1,B2_1,B3_1,B4_1,rss,ec_lda,ec_lda_rs)
524c        **** PBE96 correlation energy  corrections ****
525         t2 = t*t
526         t4 = t2*t2
527         B = -ec_lda/GAMMA
528         B = BOG/(exp(B)-1.0d0+ETA)
529         Q4 = 1.0d0 + B*t2
530         Q5 = 1.0d0 + B*t2 + B*B*t4
531         H = GAMMA*dlog(1.0d0 + BOG*Q4*t2/Q5)
532
533
534c        **** PBE96 correlation fdn and fdnc derivatives ****
535         t6   = t4*t2
536
537         B_ec = (B/BETA)*(BOG+B)
538
539         Q8  = Q5*Q5+BOG*Q4*Q5*t2
540         Q9  = 1.0d0+2*B*t2
541         H_B  = -BETA*B*t6*(2.0d0+B*t2)/Q8
542         Hrs  = H_B*B_ec*ec_lda_rs
543
544         Ht  = 2.0d0*BETA*Q9/Q8*t
545
546         ec   = ec_lda + H
547         fnc = ec  - (onethird*rs*ec_lda_rs)
548     >             - (onethird*rs*Hrs)
549     >             - (sevensixths*t*Ht)
550         fdnc = 0.5d0* Ht/ks
551
552         xce(i) = x_parameter*ex   + c_parameter*ec
553         fn(i)  = x_parameter*fnx  + c_parameter*fnc
554         fdn(i) = x_parameter*fdnx + c_parameter*fdnc
555
556
557c       write(*,*) "pbe96:",i,ec,fnc,fdnc
558
559
560      end do
561!$OMP END DO
562
563      return
564      end
565ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
566      subroutine LSDT2(a,a1,b1,b2,b3,b4,srs,ec,ec_rs)
567      real*8 a,a1,b1,b2,b3,b4,srs,ec,ec_rs
568      real*8 q0,q1,q1p,qd,ql
569      q0= -2.0d0*a*(1.0d0+a1*srs*srs)
570      q1= 2.0d0*a*srs*(b1+srs*(b2+srs*(b3+srs*b4)))
571      q1p= a*((b1/srs)+2.0d0*b2+srs*(3.0d0*b3+srs*4.0d0*b4))
572      qd=1.0d0/(q1*q1+q1)
573      ql= -dlog(qd*q1*q1)
574      ec= q0*ql
575      ec_rs= -2.0d0*a*a1*ql-q0*q1p*qd
576      return
577      end
578
579