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