1!
2! Copyright (C) 2003-2018 Quantum ESPRESSO group
3! This file is distributed under the terms of the
4! GNU General Public License. See the file `License'
5! in the root directory of the present distribution,
6! or http://www.gnu.org/copyleft/gpl.txt .
7!
8! Contributions by Eyvaz Isaev
9! Dept of Physics, Chemistry and Biology (IFM), Linkoping University, Sweden
10!
11! a) Input: Add lattice parameters units: au or Ang
12! b) Output: More info printed out
13! c) Output: Additional output file with E+PV
14!
15PROGRAM ev
16!
17!      fit of E(v) to an equation of state (EOS)
18!
19!      Interactive input:
20!         au or Ang
21!         structure
22!         equation of state
23!         input data file
24!         output data file
25!
26!      Input data file format for cubic systems:
27!         a0(1)  Etot(1)
28!         ...
29!         a0(n)  Etot(n)
30!      where a0 is the lattice parameter (a.u. or Ang)
31!      Input data file format for noncubic (e.g. hexagonal) systems:
32!         V0(1)  Etot(1)
33!         ...
34!         V0(n)  Etot(n)
35!      where V0 is the unit-cell volume (a.u.^3 or Ang^3)
36!      e.g. for an hexagonal cell,
37!         V0(i)  = sqrt(3)/2 * a^2 * c    unit-cell volume
38!         Etot(i)= min Etot(c)   for the given volume V0(i)
39!      Etot in atomic (Rydberg) units
40!
41!      Output data file format  for cubic systems:
42!      # a0=... a.u., K0=... kbar, dk0=..., d2k0=... kbar^-1, Emin=... Ry
43!      # a0=... Ang,  K0=... GPa , V0=... (a.u.)^3, V0 = Ang^3
44!         a0(1)  Etot(1) Efit(1)  Etot(1)-Efit(1)  Pfit(1)  Enth(1)
45!         ...
46!         a0(n)  Etot(n) Efit(n)  Etot(n)-Efit(n)  Pfit(n)  Enth(n)
47!      Output data file format  for noncubic systems:
48!      # V0=...(a.u.)^3, K0=... kbar, dk0=..., d2k0=... kbar^-1, Emin=... Ry
49!      # V0=...Ang^3,  K0=... GPa
50!         V0(1)  Etot(1) Efit(1)  Etot(1)-Efit(1)  Pfit(1)  Enth(1)
51!         ...
52!         V0(n)  Etot(n) Efit(n)  Etot(n)-Efit(n)  Pfit(n)  Enth(n)
53!      where
54!            a0(i), V0(i), Etot(i) as in input
55!            Efit(i) is the fitted value from the EOS
56!            Pfit(i) is the corresponding pressure from the EOS (GPa)
57!            Enth(i)=Efit(i)+Pfit(i)*V0(i) is the enthalpy (Ry)
58!!
59      USE kinds, ONLY: DP
60      USE constants, ONLY: bohr_radius_angs, ry_kbar
61      USE mp_global, ONLY : mp_startup, mp_global_end
62      USE mp_world,  ONLY : world_comm
63      USE mp,        ONLY : mp_bcast
64      USE io_global, ONLY : ionode, ionode_id
65      IMPLICIT NONE
66      INTEGER, PARAMETER:: nmaxpar=4, nmaxpt=100
67      INTEGER :: npar,npt,istat, ierr
68      CHARACTER :: bravais*8, au_unit*3, filin*256
69      REAL(DP) :: par(nmaxpar), v0(nmaxpt), etot(nmaxpt), efit(nmaxpt), &
70             fac, emin, chisq, a
71      REAL(DP), PARAMETER :: gpa_kbar = 10.0_dp
72      LOGICAL :: in_angstrom
73      CHARACTER(LEN=256) :: fileout
74  !
75  CALL mp_startup ( )
76  !
77  IF ( ionode ) THEN
78
79      WRITE(*,'(5x,"Lattice parameter or Volume are in (au, Ang) > ")', advance="NO")
80      READ(5,'(a)') au_unit
81      in_angstrom = au_unit=='Ang' .or. au_unit=='ANG' .or. &
82                    au_unit=='ang'
83      IF (in_angstrom) WRITE(*,'(5x,"Assuming Angstrom")')
84      WRITE(*,'(5x,"Enter type of bravais lattice (fcc, bcc, sc, noncubic) > ")', advance="NO")
85      READ(5,'(a)') bravais
86!
87      IF(trim(bravais)=='fcc'.or.trim(bravais)=='FCC') THEN
88         fac = 0.25d0
89      ELSEIF(trim(bravais)=='bcc'.or.trim(bravais)=='BCC') THEN
90         fac = 0.50d0
91      ELSEIF(trim(bravais)=='sc'.or.trim(bravais)=='SC') THEN
92         fac = 1.0d0
93      ELSEIF(bravais=='noncubic'.or.bravais=='NONCUBIC' .or.  &
94             trim(bravais)=='hex'.or.trim(bravais)=='HEX' ) THEN
95!         fac = sqrt(3d0)/2d0 ! not used
96         fac = 0.0_DP ! not used
97      ELSE
98         WRITE(*,'(5x,"ev: unexpected lattice <",a,">")') trim(bravais)
99         STOP
100      ENDIF
101!
102      WRITE(*,'(5x,"Enter type of equation of state :"/&
103            &  5x,"1=birch1, 2=birch2, 3=keane, 4=murnaghan > ")', advance="NO")
104      READ(5,*) istat
105      IF(istat==1 .or. istat==4) THEN
106         npar=3
107      ELSEIF(istat==2 .or. istat==3) THEN
108         npar=4
109      ELSE
110         WRITE(*,'(5x,"Unexpected eq. of state ",i2)') istat
111         STOP
112      ENDIF
113      WRITE(*,'(5x,"Input file > ")', advance="NO")
114      READ(5,'(a)') filin
115      OPEN(unit=2,file=filin,status='old',form='formatted',iostat=ierr)
116      IF (ierr/=0) THEN
117         ierr= 1
118         GO TO 99
119      END IF
120  10  CONTINUE
121      emin=1d10
122      DO npt=1,nmaxpt
123         IF (bravais=='noncubic'.or.bravais=='NONCUBIC' .or. &
124             bravais=='hex'.or.bravais=='HEX' ) THEN
125            READ(2,*,err=10,END=20) v0(npt), etot(npt)
126            IF (in_angstrom) v0(npt)=v0(npt)/bohr_radius_angs**3
127         ELSE
128            READ(2,*,err=10,END=20) a, etot(npt)
129            IF (in_angstrom) a = a/bohr_radius_angs
130            v0  (npt) = fac*a**3
131         ENDIF
132         IF(etot(npt)<emin) THEN
133            par(1) = v0(npt)
134            emin = etot(npt)
135         ENDIF
136      ENDDO
137
138      npt = nmaxpt+1
139  20  npt = npt-1
140!
141! par(1) = V, Volume of the unit cell in (a.u.^3)
142! par(2) = B, Bulk Modulus (in KBar)
143! par(3) = dB/dP (adimensional)
144! par(4) = d^2B/dP^2 (in KBar^(-1), used only by 2nd order formulae)
145!
146      par(2) = 500.0d0
147      par(3) = 5.0d0
148      par(4) = -0.01d0
149
150      CALL find_minimum &
151           (npar,par,chisq)
152!
153      CALL write_results &
154           (npt,in_angstrom,fac,v0,etot,efit,istat,par,npar,emin,chisq, &
155            fileout)
156!
157      CALL write_evdata_xml  &
158           (npt,fac,v0,etot,efit,istat,par,npar,emin,chisq,fileout, ierr)
159
160      IF (ierr /= 0) GO TO 99
161    ENDIF
16299  CALL mp_bcast ( ierr, ionode_id, world_comm )
163    IF ( ierr == 1) THEN
164       CALL errore( 'ev', 'file '//trim(filin)//' cannot be opened', ierr )
165    ELSE IF ( ierr == 2 ) THEN
166       CALL errore( 'ev', 'file '//trim(fileout)//' cannot be opened', ierr )
167    ELSE IF ( ierr == 11 ) THEN
168       CALL errore( 'write_evdata_xml', 'no free units to write ', ierr )
169    ELSE IF ( ierr == 12 ) THEN
170       CALL errore( 'write_evdata_xml', 'error opening the xml file ', ierr )
171    ENDIF
172
173  CALL mp_global_end()
174
175      STOP
176    CONTAINS
177!
178!-----------------------------------------------------------------------
179      SUBROUTINE eqstate(npar,par,chisq)
180!-----------------------------------------------------------------------
181!
182      IMPLICIT NONE
183      INTEGER, INTENT(in) :: npar
184      REAL(DP), INTENT(in) :: par(npar)
185      REAL(DP), INTENT(out):: chisq
186      INTEGER :: i
187      REAL(DP) :: k0, dk0, d2k0, c0, c1, x, vol0, ddk
188!
189      vol0 = par(1)
190      k0   = par(2)/ry_kbar ! converts k0 to Ry atomic units...
191      dk0  = par(3)
192      d2k0 = par(4)*ry_kbar ! and d2k0/dp2 to (Ry a.u.)^(-1)
193!
194      IF(istat==1.or.istat==2) THEN
195         IF(istat==1) THEN
196            c0 = 0.0d0
197         ELSE
198            c0 = ( 9.d0*k0*d2k0 + 9.d0*dk0**2-63.d0*dk0+143.d0 )/48.d0
199         ENDIF
200         c1 = 3.d0*(dk0-4.d0)/8.d0
201         DO i=1,npt
202            x = vol0/v0(i)
203            efit(i) = 9.d0*k0*vol0*( (-0.5d0+c1-c0)*x**(2.d0/3.d0)/2.d0 &
204                         +( 0.50-2.d0*c1+3.d0*c0)*x**(4.d0/3.d0)/4.d0 &
205                         +(       c1-3.d0*c0)*x**(6.d0/3.d0)/6.d0 &
206                         +(            c0)*x**(8.d0/3.d0)/8.d0 &
207                         -(-1.d0/8.d0+c1/6.d0-c0/8.d0) )
208         ENDDO
209      ELSE
210         IF(istat==3) THEN
211            ddk = dk0 + k0*d2k0/dk0
212         ELSE
213            ddk = dk0
214         ENDIF
215         DO i=1,npt
216            efit(i) = - k0*dk0/ddk*vol0/(ddk-1.d0) &
217            + v0(i)*k0*dk0/ddk**2*( (vol0/v0(i))**ddk/(ddk-1.d0)+1.d0) &
218            - k0*(dk0-ddk)/ddk*( v0(i)*log(vol0/v0(i)) + v0(i)-vol0 )
219         ENDDO
220      ENDIF
221!
222!      emin = equilibrium energy obtained by minimizing chi**2
223!
224      emin = 0.0d0
225      DO i = 1,npt
226         emin = emin + etot(i)-efit(i)
227      ENDDO
228      emin = emin/npt
229!
230      chisq = 0.0d0
231      DO i = 1,npt
232          efit(i) = efit(i)+emin
233          chisq   = chisq + (etot(i)-efit(i))**2
234      ENDDO
235      chisq = chisq/npt
236!
237      RETURN
238    END SUBROUTINE eqstate
239!
240!-----------------------------------------------------------------------
241      SUBROUTINE write_results &
242            (npt,in_angstrom,fac,v0,etot,efit,istat,par,npar,emin,chisq, &
243             filout)
244!-----------------------------------------------------------------------
245!
246      IMPLICIT NONE
247      INTEGER, INTENT(in) :: npt, istat, npar
248      REAL(DP), INTENT(in):: v0(npt), etot(npt), efit(npt), emin, chisq, fac
249      REAL(DP), INTENT(inout):: par(npar)
250      REAL(DP), EXTERNAL :: keane, birch
251      LOGICAL, INTENT(in) :: in_angstrom
252      CHARACTER(len=256), intent(inout) :: filout
253      !
254      REAL(DP) :: p(npt), epv(npt)
255      INTEGER :: i, iun
256      LOGICAL :: exst
257
258      WRITE(*,'(5x,"Output file > ")', advance="NO")
259      READ (5,'(a)') filout
260      IF(filout/=' ') THEN
261         iun=8
262         INQUIRE(file=filout,exist=exst)
263         IF (exst) WRITE(*,'(5x,"Beware: file ",A," will be overwritten")')&
264                  trim(filout)
265         OPEN(unit=iun,file=filout,form='formatted',status='unknown', &
266              iostat=ierr)
267         IF (ierr/=0) THEN
268            ierr= 2
269            GO TO 99
270         END IF
271      ELSE
272         iun=6
273      ENDIF
274
275      IF(istat==1) THEN
276         WRITE(iun,'("# equation of state: birch 1st order.  chisq = ", &
277                   & d12.4)') chisq
278      ELSEIF(istat==2) THEN
279         WRITE(iun,'("# equation of state: birch 3rd order.  chisq = ", &
280                   & d12.4)') chisq
281      ELSEIF(istat==3) THEN
282         WRITE(iun,'("# equation of state: keane.            chisq = ", &
283                   & d12.4)') chisq
284      ELSEIF(istat==4) THEN
285         WRITE(iun,'("# equation of state: murnaghan.        chisq = ", &
286                   & d12.4)') chisq
287      ENDIF
288
289      IF(istat==1 .or. istat==4) par(4) = 0.0d0
290
291      IF(istat==1 .or. istat==2) THEN
292         DO i=1,npt
293            p(i)=birch(v0(i)/par(1),par(2),par(3),par(4))
294         ENDDO
295      ELSE
296         DO i=1,npt
297            p(i)=keane(v0(i)/par(1),par(2),par(3),par(4))
298         ENDDO
299      ENDIF
300
301      DO i=1,npt
302         epv(i) = etot(i) + p(i)*v0(i) / ry_kbar
303      ENDDO
304
305      IF ( fac /= 0.0_dp ) THEN
306! cubic case
307         WRITE(iun,'("# a0 =",f8.4," a.u., k0 =",i5," kbar, dk0 =", &
308                    &f6.2," d2k0 =",f7.3," emin =",f11.5)') &
309            (par(1)/fac)**(1d0/3d0), int(par(2)), par(3), par(4), emin
310         WRITE(iun,'("# a0 =",f9.5," Ang, k0 =", f6.1," GPa,  V0 = ", &
311                  & f8.2," (a.u.)^3,  V0 =", f8.2," A^3 ",/)') &
312           & (par(1)/fac)**(1d0/3d0)*bohr_radius_angs, par(2)/gpa_kbar, &
313             par(1), par(1)*bohr_radius_angs**3
314
315        WRITE(iun,'(73("#"))')
316        WRITE(iun,'("# Lat.Par", 7x, "E_calc", 8x, "E_fit", 7x, &
317             & "E_diff", 4x, "Pressure", 6x, "Enthalpy")')
318        IF (in_angstrom) THEN
319           WRITE(iun,'("# Ang", 13x, "Ry", 11x, "Ry", 12x, &
320             & "Ry", 8x, "GPa", 11x, "Ry")')
321           WRITE(iun,'(73("#"))')
322           WRITE(iun,'(f9.5,2x,f12.5, 2x,f12.5, f12.5, 3x, f8.2, 3x,f12.5)') &
323              & ( (v0(i)/fac)**(1d0/3d0)*bohr_radius_angs, etot(i), efit(i),  &
324              & etot(i)-efit(i), p(i)/gpa_kbar, epv(i), i=1,npt )
325        ELSE
326           WRITE(iun,'("# a.u.",12x, "Ry", 11x, "Ry", 12x, &
327             & "Ry", 8x, "GPa", 11x, "Ry")')
328           WRITE(iun,'(73("#"))')
329           WRITE(iun,'(f9.5,2x,f12.5, 2x,f12.5, f12.5, 3x, f8.2, 3x,f12.5)') &
330              & ( (v0(i)/fac)**(1d0/3d0), etot(i), efit(i),  &
331              & etot(i)-efit(i), p(i)/gpa_kbar, epv(i), i=1,npt )
332        ENDIF
333
334      ELSE
335! noncubic case
336         WRITE(iun,'("# V0 =",f8.2," a.u.^3,  k0 =",i5," kbar,  dk0 =", &
337                    & f6.2,"  d2k0 =",f7.3,"  emin =",f11.5)') &
338                    & par(1), int(par(2)), par(3), par(4), emin
339
340         WRITE(iun,'("# V0 =",f8.2,"  Ang^3,  k0 =",f6.1," GPa"/)') &
341                    & par(1)*bohr_radius_angs**3, par(2)/gpa_kbar
342
343        WRITE(iun,'(74("#"))')
344        WRITE(iun,'("# Vol.", 8x, "E_calc", 8x, "E_fit", 7x, &
345             & "E_diff", 4x, "Pressure", 6x, "Enthalpy")')
346        IF (in_angstrom) THEN
347          WRITE(iun,'("# Ang^3", 9x, "Ry", 11x, "Ry", 12x, &
348             & "Ry", 8x, "GPa", 11x, "Ry")')
349          WRITE(iun,'(74("#"))')
350           WRITE(iun,'(f8.2,2x,f12.5, 2x,f12.5, f12.5, 3x, f8.2, 3x,f12.5)') &
351              ( v0(i)*bohr_radius_angs**3, etot(i), efit(i),  &
352               etot(i)-efit(i), p(i)/gpa_kbar, epv(i), i=1,npt )
353         else
354          WRITE(iun,'("# a.u.^3",8x, "Ry", 11x, "Ry", 12x, &
355             & "Ry", 8x, "GPa", 11x, "Ry")')
356          WRITE(iun,'(74("#"))')
357           WRITE(iun,'(f8.2,2x,f12.5, 2x,f12.5, f12.5, 3x, f8.2, 3x,f12.5)') &
358              ( v0(i), etot(i), efit(i),  &
359               etot(i)-efit(i), p(i)/gpa_kbar, epv(i), i=1,npt )
360         end if
361
362      ENDIF
363      IF(filout/=' ') CLOSE(unit=iun)
364 99   RETURN
365    END SUBROUTINE write_results
366!
367
368      ! This subroutine is passed to LMDIF to be minimized
369      SUBROUTINE SCHISQ(m_, n_, par_, f_, i_)
370           IMPLICIT NONE
371           INTEGER,INTENT(in)  :: m_, n_
372           INTEGER,INTENT(inout)   :: i_
373           REAL(DP),INTENT(in)    :: par_(n_)
374           REAL(DP),INTENT(out)   :: f_(m_)
375           REAL(DP) :: chisq_
376           IF(m_/=n_) THEN
377              i_ = -999
378              RETURN
379           ENDIF
380           !
381           CALL eqstate(n_,par_,chisq_)
382           f_=0._dp
383           f_(1) = chisq_
384        END SUBROUTINE
385
386!-----------------------------------------------------------------------
387      SUBROUTINE find_minimum(npar,par,chisq)
388!-----------------------------------------------------------------------
389!
390      USE lmdif_module, ONLY : lmdif1
391      IMPLICIT NONE
392      INTEGER ,INTENT(in)  :: npar
393      REAL(DP),INTENT(out) :: par(nmaxpar)
394      REAL(DP),INTENT(out) :: chisq
395      !
396      REAL(DP) :: xi(npar,npar), vchisq(npar)
397      INTEGER :: i
398      INTEGER :: lwa, iwa(npar)
399      REAL(DP),ALLOCATABLE :: wa(:)
400      !
401      xi = 0._dp
402      FORALL(i=1:npar) xi(i,i) = 1._dp
403      par(1) = v0(npt/2)
404      par(2) = 500.0d0
405      par(3) = 5.0d0
406      par(4) = -0.01d0 ! unused for some eos
407
408      lwa = npar**2 + 6*npar
409      ALLOCATE(wa(lwa))
410      !
411      CALL lmdif1(SCHISQ, npar, npar, par, vchisq, 1.d-12, i, iwa, wa, lwa)
412      DEALLOCATE(wa)
413      !
414      IF(i>0 .and. i<5) THEN
415         PRINT*, "Minimization succeeded"
416      ELSEIF(i>=5) THEN
417         PRINT*, "Minimization stopped before convergence"
418      ELSEIF(i<=0) THEN
419        PRINT*, "Minimization error"
420        STOP
421      ENDIF
422      !chisq = vchisq(1)
423      !
424      CALL eqstate(npar,par,chisq)
425
426      END SUBROUTINE
427!-----------------------------------------------------------------------
428
429  END PROGRAM ev
430
431      FUNCTION birch(x,k0,dk0,d2k0)
432!
433      USE kinds, ONLY : DP
434      IMPLICIT NONE
435      REAL(DP) birch, x, k0,dk0, d2k0
436      REAL(DP) c0, c1
437
438      IF(d2k0/=0.d0) THEN
439         c0 = (9.d0*k0*d2k0 + 9.d0*dk0**2 - 63.d0*dk0 + 143.d0 )/48.d0
440      ELSE
441         c0 = 0.0d0
442      ENDIF
443      c1 = 3.d0*(dk0-4.d0)/8.d0
444      birch = 3.d0*k0*( (-0.5d0+  c1-  c0)*x**( -5.d0/3.d0) &
445           +( 0.5d0-2.d0*c1+3.0d0*c0)*x**( -7.d0/3.d0) &
446           +(       c1-3*c0)*x**( -9.0d0/3d0) &
447           +(            c0)*x**(-11.0d0/3d0) )
448      RETURN
449    END FUNCTION birch
450!
451      FUNCTION keane(x,k0,dk0,d2k0)
452!
453      USE kinds, ONLY : DP
454      IMPLICIT NONE
455      REAL(DP) keane, x, k0, dk0, d2k0, ddk
456
457      ddk = dk0 + k0*d2k0/dk0
458      keane = k0*dk0/ddk**2*( x**(-ddk) - 1d0 ) + (dk0-ddk)/ddk*log(x)
459
460      RETURN
461    END FUNCTION keane
462!-----------------------------------------------------------------------
463    SUBROUTINE write_evdata_xml &
464        (npt,fac,v0,etot,efit,istat,par,npar,emin,chisq,filout, ierr)
465!-----------------------------------------------------------------------
466!
467  USE kinds,     ONLY : dp
468  USE constants, ONLY : ry_kbar, bohr_radius_angs
469  IMPLICIT NONE
470  INTEGER, INTENT(in) :: npt, istat, npar
471  REAL(DP), INTENT(in):: v0(npt), etot(npt), efit(npt), emin, chisq, fac
472  REAL(DP), INTENT(in):: par(npar)
473  CHARACTER(len=256), INTENT(IN) :: filout
474  INTEGER, INTENT(out) :: ierr
475  !
476  INTEGER :: iunout = 11
477  REAL(DP) :: p(npt), volume(2), a0(2), alldata(6,npt)
478  INTEGER :: i, iun
479  CHARACTER(len=256) :: filename
480  REAL(DP), EXTERNAL :: birch, keane
481
482  IF (filout/=' ') THEN
483     filename = TRIM(filout) // '.xml'
484  ELSE
485     filename = 'ev.xml'
486  ENDIF
487  !
488  ! ... open XML descriptor
489  !
490  OPEN ( UNIT=iunout, FILE = TRIM( filename ), FORM='formatted', IOSTAT = ierr )
491  IF ( ierr /= 0 ) THEN
492     WRITE (6,*) 'Failed opening file ' // TRIM(filename)
493     RETURN
494  END IF
495
496  WRITE (iunout,'("<xml>")')
497  WRITE (iunout,'("<EQUATIONS_OF_STATE>")')
498  WRITE (iunout,'("<EQUATION_TYPE>")')
499  IF (istat==1) THEN
500     WRITE (iunout,'("Birch 1st order")')
501  ELSEIF (istat==2) THEN
502     WRITE (iunout,'("Birch 2nd order")')
503  ELSEIF (istat==3) THEN
504     WRITE (iunout,'("Keane")')
505  ELSEIF (istat==4) THEN
506     WRITE (iunout,'("Murnaghan")')
507  ENDIF
508  WRITE (iunout,'("</EQUATION_TYPE>")')
509  WRITE (iunout,'("<CHI_SQUARE>")')
510  WRITE (iunout,'(1pe25.12)') chisq
511  WRITE (iunout,'("</CHI_SQUARE>")')
512  WRITE (iunout,'("</EQUATIONS_OF_STATE>")')
513
514  IF (istat==1 .or. istat==2) THEN
515     DO i=1,npt
516        p(i)=birch(v0(i)/par(1),par(2),par(3),par(4))
517     ENDDO
518  ELSE
519     DO i=1,npt
520        p(i)=keane(v0(i)/par(1),par(2),par(3),par(4))
521     ENDDO
522  ENDIF
523
524  DO i=1,npt
525     alldata (1,i) = v0(i)
526     alldata (2,i) = etot(i)
527     alldata (3,i) = efit(i)
528     alldata (4,i) = etot(i) - efit(i)
529     alldata (5,i) = p(i)
530     alldata (6,i) = etot(i) + p(i) * v0(i) / ry_kbar
531  ENDDO
532
533  WRITE (iunout,'("<EQUATIONS_PARAMETERS>")')
534
535  volume(1)=par(1)
536  volume(2)=par(1)*bohr_radius_angs**3
537  WRITE (iunout, '("<EQUILIBRIUM_VOLUME_AU_A>")')
538  WRITE (iunout, '(2(1pe25.15))') volume(:)
539  WRITE (iunout, '("</EQUILIBRIUM_VOLUME_AU_A>")')
540  WRITE (iunout, '("<BULK_MODULUS_KBAR>")')
541  WRITE (iunout, '(1pe25.15)') par(2)
542  WRITE (iunout, '("</BULK_MODULUS_KBAR>")')
543  WRITE (iunout, '("<DERIVATIVE_BULK_MODULUS>")')
544  WRITE (iunout, '(1pe25.15)') par(3)
545  WRITE (iunout, '("</DERIVATIVE_BULK_MODULUS>")')
546  WRITE (iunout, '("<SECOND_DERIVATIVE_BULK_MODULUS>")')
547  WRITE (iunout, '(1pe25.15)') par(4)
548  WRITE (iunout, '("</SECOND_DERIVATIVE_BULK_MODULUS>")')
549  WRITE (iunout, '("<MINIMUM_ENERGY_RY>")')
550  WRITE (iunout, '(1pe25.15)') emin
551  WRITE (iunout, '("</MINIMUM_ENERGY_RY>")')
552  WRITE (iunout, '("<CELL_FACTOR>")')
553  WRITE (iunout, '(1pe25.15)') fac
554  WRITE (iunout, '("</CELL_FACTOR>")')
555  IF (fac /= 0.0_DP) THEN
556     a0(1) = (par(1)/fac)**(1d0/3d0)
557     a0(2) = (par(1)/fac)**(1d0/3d0) * bohr_radius_angs
558     WRITE (iunout, '("<CELL_PARAMETER_AU_A>")')
559     WRITE (iunout, '(2(1pe25.15))') a0
560     WRITE (iunout, '("</CELL_PARAMETER_AU_A>")')
561  ENDIF
562  WRITE (iunout,'("</EQUATIONS_PARAMETERS>")')
563
564  WRITE (iunout,'("<FIT_CHECK>")')
565  WRITE (iunout,'("<NUMBER_OF_DATA>")')
566  WRITE (iunout,'(i8)') npt
567  WRITE (iunout,'("</NUMBER_OF_DATA>")')
568  WRITE (iunout,'("<VOL_ENE_EFIT_DELTA_P_GIBBS>")')
569  WRITE (iunout,'(6(1pe25.15))') alldata(:,:)
570  WRITE (iunout,'("</VOL_ENE_EFIT_DELTA_P_GIBBS>")')
571
572  WRITE (iunout,'("</FIT_CHECK>")')
573  CLOSE (unit=iunout, status='keep')
574
575  RETURN
576END SUBROUTINE write_evdata_xml
577
578