1! --- 2! Copyright (C) 1996-2016 The SIESTA group 3! This file is distributed under the terms of the 4! GNU General Public License: see COPYING in the top directory 5! or http://www.gnu.org/copyleft/gpl.txt . 6! See Docs/Contributors.txt for a list of contributors. 7! --- 8subroutine read_options( na, ns, nspin ) 9 10 ! Subroutine to read the options for the SIESTA program 11 ! 12 ! It uses the FDF (Flexible Data Format) package 13 ! of J.M.Soler and A.Garcia 14 ! 15 ! Writen by P.Ordejon, December'96 16 ! Modified for introduction of dynamic memory in SIESTA by JDG Sept 99 17 ! Wrapping of most fdf and broadcast calls: A. Garcia, June 2005 18 ! 19 20 use siesta_options 21 use precision, only : dp, grid_p 22 use parallel, only : IOnode, Nodes 23 use fdf 24 use files, only : slabel 25 use files, only : filesOut_t ! derived type for output file names 26 use sys 27 use units, only : eV, Ang, Kelvin 28 use siesta_cml 29 use m_target_stress, only: set_target_stress 30 use m_spin, only: print_spin_options 31 32 use m_charge_add, only : read_charge_add 33 use m_hartree_add, only : read_hartree_add 34 35 use m_mixing_scf, only: mixers_scf_init 36 use m_mixing_scf, only: mixers_scf_print, mixers_scf_print_block 37 38 use m_cite, only: add_citation 39 40 implicit none 41 !----------------------------------------------------------- Input Variables 42 ! integer na : Number of atoms 43 ! integer ns : Number of species 44 ! integer nspin : Number of spin-components 45 ! 1=non-polarized, 2=polarized, 4=non-collinear, 46 ! 8=spin-orbit 47 48 integer, intent(in) :: na, ns, nspin 49 50 ! This routine sets variables in the 'siesta_options' module 51 52 ! The following are comment lines that should be merged into 'siesta_options'. 53 54 ! real*8 charnet : Net charge (in units of |e|) 55 ! logical outlng : Long (true) or short (false) output 56 ! real*8 g2cut : PW cutoff energy (Ry) 57 ! logical negl : True = Neglect interactions between 58 ! non-overlaping orbitals (coming from 59 ! KB projectors) 60 ! integer nscf : Maximum number of SCF cycles per time step 61 ! real*8 dDtol : Maximum Density Matrix tolerance in SCF 62 ! real*8 Energy_tolerance : Maximum Total energy tolerance in SCF 63 ! real*8 Harris_tolerance : Maximum Harris energy tolerance in SCF 64 ! logical mix : Perform mix in first SCF step 65 ! real*8 wmix : Amount of output DM for new DM 66 ! integer isolve : Method of solution. 0 = Diagonalization 67 ! 1 = Order-N 68 ! 2 = Transiesta 69 ! 3 = OMM 70 ! 4 = PEXSI 71 ! 5 = (Matrix write) 72 ! real*8 temp : Temperature for Fermi smearing (Ry) 73 ! logical fixspin : Fix the spin of the system? 74 ! real*8 total_spin : Total spin of the system 75 ! integer ncgmax : Maximum number of CG steps for 76 ! band structure energy minimization 77 ! real*8 etol : Relative tolerance in CG minimization 78 ! of band structure energy 79 ! real*8 eta(2) : Fermi level parameter of Kim functional 80 ! real*8 rcoor : Cutoff radius of LWF's (Bohr) 81 ! integer ioptlwf : Option to build LWF's according to: 82 ! 0 = Read blindly from disk 83 ! 1 = Functional of Kim et al. 84 ! 2 = Functional of Ordejon-Mauri 85 ! logical chebef : Compute the chemical potential 86 ! logical noeta : Use computed Chem.pot. instead of eta 87 ! real*8 rcoorcp : Cutoff (Bohr) to compute the chem.pot. 88 ! real*8 beta : Inverse temperature to compute chem.pot. 89 ! integer pmax : Order of Chebi expansion for chem.pot. 90 ! integer idyn : Atomic dynamics option: 91 ! 0 = Geometry optimization 92 ! 1 = Standard MD run (Verlet) 93 ! 2 = Nose thermostat MD 94 ! 3 = Parrinello-Rahman MD 95 ! 4 = Nose thermostat + Parrinello-Rahman MD 96 ! 5 = Annealing MD 97 ! 6 = Force constants 98 ! 7 = Deprecated (Forces for PHONON program) 99 ! 8 = Force evaluation 100 ! 9 = Explicit set of coordinates 101 ! 10 = Lua controlled dynamics 102 ! integer istart : Initial time step for MD 103 ! integer ifinal : Final time step for MD 104 ! integer nmove : Number of steps in *any* MD/optimization 105 ! real*8 ftol : Maximum force for structural optimization 106 ! real*8 strtol : Maximum stress for structural optimization 107 ! integer ianneal : Annealing option for idyn = 5 108 ! 1 = Temperature 109 ! 2 = Pressure 110 ! 3 = Temperature and Pressure 111 ! integer iquench : Quench option: 0 = No; 1 = Yes; 2 = Fire 112 ! real*8 dt : Length of time step (fs) 113 ! real*8 dx : Atomic displacement for Force Constants 114 ! calculation 115 ! integer ia1 : First atom to displace for force constants 116 ! integer ia2 : Last atom to displace for force constants 117 ! real*8 dxmax : Maximum atomic displacement in one atomic move 118 ! real*8 tt : Target temperature (Kelvin) 119 ! real*8 tp : Target Pressure (Ry/Bohr**3) 120 ! real*8 mn : Mass of Nose variable (Ry/fs**2) 121 ! real*8 mpr : Mass of Parrinello-R. variable (Ry/fs**2) 122 ! real*8 bulkm : Estimate of bulk modulus (Ry/Bohr**3) 123 ! real*8 taurelax : Annealing time to reach targer T and P (fs) 124 ! logical usesavelwf : True = try to use continuation LWF files 125 ! from disk 126 ! logical usesavedm : True = try to use continuation DM files 127 ! from disk 128 ! logical usesavecg : True = try to use continuation CG files 129 ! from disk 130 ! integer mullipop : Option for Mulliken Pop. analysis 131 ! logical init_anti_ferro : Spin initialization for spin-polarized 132 ! .true. -> Antiferro 133 ! .false. -> Ferro 134 ! integer maxsav : Number of density-matrices stored for Pulay 135 ! mixing. .lt.2 => linear mixing only 136 ! .ge.2 => pulay mixing 137 ! integer nkick : Perform a linear mixing eack nkick scf cycles 138 ! real*8 wmixkick : Mixing parameter for linear mixing each nkick scf 139 ! cycles 140 ! logical pulfile : Use file (.true.) or memory (.false.) 141 ! to store Pulay miximg intermediate vectors 142 ! Default: .false. 143 ! real*8 tempinit : Initial temperature (Kelvin) of the MD simulation 144 ! logical dumpcharge : True: Dump information to plot charge contours 145 ! by the external DENCHAR application program. 146 ! (This is now obsolete: info will appear in .RHO file) 147 ! logical varcel : variable shape for optimization or dynamics 148 ! logical harrisfun : swith that indicates if harris functional will 149 ! be used or not 150 ! real*8 occtol : Occupancy threshold for DM build 151 ! integer broyden_maxit : Number of histories saved in Broyden SCF mixing 152 ! logical require_energy_convergence : Impose E. conv. criterion? 153 ! logical broyden_optim : Use Broyden method for optimization 154 ! logical want_domain_decomposition: Use domain decomposition for orbitals in O(N) 155 ! logical want_spatial_decomposition: Use spatial decomposition for orbitals in O(N) 156 157 158 !----------------------------------------------------------- Local Variables 159 real(dp) :: tcp 160 161 character(len=22) :: annop, dyntyp 162 character(len=13) :: lwfopt 163 character(len=30) :: ctmp 164 character(len=6) :: method 165 166 logical :: qnch, qnch2 167 logical :: tBool 168 169 !--------------------------------------------------------------------- BEGIN 170 ! New template, using fdf 171 ! 172 ! param = fdf_get('ParamName', param_default) 173 ! if (ionode) write(6,'(a,i)'), 174 ! . 'redata: ParamName = ',param 175 ! if (cml_p) call cmlAddParameter(xf=mainXML, name='ParamName', 176 ! . value=param, dictref='siesta:param') 177 178 ! 179 ! cml_p is only true in the master node 180 ! 181 ! Start of code 182 183 if (cml_p) then 184 call cmlStartParameterList(mainXML, title='Input Parameters') 185 endif 186 187 ! for cml output, find the system name & label 188 if (cml_p) then 189 call cmlAddParameter(xf=mainXML, name='SystemName', & 190 value=trim(sname), dictref='siesta:sname') 191 call cmlAddParameter(xf=mainXML, name='SystemLabel', & 192 value=trim(slabel), dictref='siesta:slabel') 193 endif 194 195 ! Start by printing out spin-configuration 196 call print_spin_options() 197 198 ! H setup only 199 h_setup_only = fdf_get('HSetupOnly', .false.) 200 if (ionode .and. h_setup_only) then 201 write(6,1) 'redata: H Setup Only', h_setup_only 202 endif 203 204 ! Type of output 205 outlng = fdf_get('LongOutput', .false.) 206 if (ionode) then 207 write(6,1) 'redata: Long output', outlng 208 endif 209 210 if (cml_p) then 211 call cmlAddParameter( xf=mainXML, name='LongOutput', & 212 value=outlng, dictRef='siesta:verbosity' ) 213 endif 214 215 ! Write about Number of species, as before 216 if (ionode) then 217 write(6,4) 'redata: Number of Atomic Species', ns 218 endif 219 220 if (ns .le. 0) then 221 call die( 'redata: ERROR: Number of species must be larger than zero.' ) 222 endif 223 224 if (cml_p) then 225 call cmlAddParameter( xf=mainXML, title='NumberOfSpecies', & 226 value=ns, dictRef='siesta:ns', units="cmlUnits:countable" ) 227 endif 228 229 ! Dump information to plot charge contours 230 ! by the external DENCHAR application program. 231 dumpcharge = fdf_get('WriteDenchar',.false.) 232 if (ionode) then 233 write(6,2) 'redata: Charge density info will appear in .RHO file' 234 endif 235 236 if (cml_p) then 237 call cmlAddParameter( xf=mainXML, name='WriteDenChar', value=dumpcharge) 238 endif 239 240 ! Perform Mulliken Population Analysis 241 mullipop = fdf_get('WriteMullikenPop', 0) 242 if (mullipop == 0 .and. outlng) then 243 mullipop = 1 244 endif 245 ! <L> output 246 orbmoms = fdf_get( 'WriteOrbMom' , .false. ) 247 248 if (ionode) then 249 select case (mullipop) 250 case(0) 251 write(6,3) 'redata: Write Mulliken Pop.','NO' 252 case(1) 253 write(6,3) 'redata: Write Mulliken Pop.','Atomic and Orbital charges' 254 case(2) 255 write(6,3)'redata: Write Mulliken Pop.','Atomic and Orbital charges' 256 write(6,10) 'plus Atomic Overlap Populations' 257 case(3) 258 write(6,3)'redata: Write Mulliken Pop.','Atomic and Orbital charges' 259 write(6,10) 'plus Atomic Overlap Populations' 260 write(6,10) 'plus Orbital Overlap Populations' 261 case default 262 call die( 'redata: Invalid value for WriteMullikenPop' ) 263 end select 264 endif 265 266 if (cml_p) then 267 call cmlAddParameter( xf=mainXML, name='WriteMullikenPop', value=mullipop, & 268 units="cmlUnits:dimensionless" ) 269 endif 270 271 ! Perform Hirshfeld and/or Voronoi Population Analysis 272 hirshpop= fdf_get('WriteHirshfeldPop',.false.) 273 voropop= fdf_get('WriteVoronoiPop',.false.) 274 partial_charges_at_every_geometry = & 275 fdf_get('PartialChargesAtEveryGeometry',.false.) 276 partial_charges_at_every_scf_step = & 277 fdf_get('PartialChargesAtEveryScfStep',.false.) 278 279 280 if ( fdf_get('Compat.Matel.NRTAB', .false.) ) then 281 matel_NRTAB = 128 282 else 283 matel_NRTAB = 1024 284 end if 285 if ( IONode ) then 286 write(6,4) 'redata: Matel table size (NRTAB)', matel_NRTAB 287 end if 288 if (cml_p) then 289 call cmlAddParameter( xf=mainXML, name='MatelNRTAB',value=matel_NRTAB, & 290 dictRef='siesta:matel_nrtab', units="cmlUnits:countable") 291 end if 292 293 ! Planewave cutoff of the real space mesh ... 294 g2cut = fdf_get('MeshCutoff',300._dp,'Ry') 295 if (ionode) then 296 write(6,6) 'redata: Mesh Cutoff', g2cut,' Ry' 297 endif 298 299 if (cml_p) then 300 call cmlAddParameter( xf=mainXML, name='MeshCutOff', value=g2cut, & 301 dictRef='siesta:g2max', units='siestaUnits:Ry' ) 302 endif 303 304 ! Net charge in the cell ... 305 charnet = fdf_get('NetCharge',0.0_dp) 306 if (ionode) then 307 write(6,6) 'redata: Net charge of the system',charnet,' |e|' 308 endif 309 310 if (cml_p) then 311 call cmlAddParameter( xf=mainXML, name='NetCharge', value=charnet, & 312 dictRef='siesta:NetCharge', units='siestaUnits:e__') 313 endif 314 315 ! SCF Loop parameters ... 316 ! Minimum/Maximum number of SCF iterations 317 min_nscf = fdf_get('MinSCFIterations',0) 318 nscf = fdf_get('MaxSCFIterations',1000) 319 SCFMustConverge = fdf_get('SCF.MustConverge', .true.) 320 if (ionode) then 321 write(6,4) 'redata: Min. number of SCF Iter',min_nscf 322 write(6,4) 'redata: Max. number of SCF Iter',nscf 323 if (SCFMustConverge) then 324 write(6,2) 'redata: SCF convergence failure will abort job' 325 endif 326 endif 327 328 if (cml_p) then 329 call cmlAddParameter( xf=mainXML, name='MaxSCFIterations', & 330 value=nscf, dictRef='siesta:maxscf', & 331 units="cmlUnits:countable") 332 call cmlAddParameter( xf=mainXML, name='MinSCFIterations', & 333 value=min_nscf, dictRef='siesta:minscf', & 334 units="cmlUnits:countable") 335 endif 336 337 call fdf_deprecated('TS.MixH','SCF.Mix') 338 call fdf_deprecated('MixHamiltonian','SCF.Mix') 339 call fdf_deprecated('MixCharge','SCF.Mix') 340 341 ! Note, since 4.1 mixing the Hamiltonian is the default option! 342 mixH = fdf_get('TS.MixH',.true.) ! Catch old-style keyword (prefer new key) 343 mixH = fdf_get('MixHamiltonian',mixH) 344 mix_charge = fdf_get('MixCharge',.false.) 345 346 if ( mix_charge ) then 347 ctmp = 'charge' 348 else if ( mixH ) then 349 ctmp = 'Hamiltonian' 350 else 351 ctmp = 'density' 352 end if 353 354 ctmp = fdf_get('SCF.Mix', trim(ctmp)) 355 if ( leqi(ctmp, 'charge') .or. & 356 leqi(ctmp, 'rho') ) then 357 mix_charge = .true. 358 mixH = .false. 359 else if ( leqi(ctmp, 'Hamiltonian') & 360 .or. leqi(ctmp, 'H') ) then 361 mix_charge = .false. 362 mixH = .true. 363 else if ( leqi(ctmp, 'density') & 364 .or. leqi(ctmp, 'density-matrix') & 365 .or. leqi(ctmp, 'DM') ) then 366 mix_charge = .false. 367 mixH = .false. 368 else 369 call die('Unrecognized option for: SCF.Mix. Please see the manual.') 370 end if 371 372 if ( IONode ) then 373 if ( mix_charge ) then 374 write(6,3) 'redata: SCF mix quantity', 'charge' 375 else if ( mixH ) then 376 write(6,3) 'redata: SCF mix quantity', 'Hamiltonian' 377 else 378 write(6,3) 'redata: SCF mix quantity', 'density-matrix' 379 end if 380 end if 381 382 ! Options for pre-4.0 compatibility 383 compat_pre_v4_DM_H = fdf_get('Compat-pre-v4-DM-H',.false.) 384 mix_after_convergence = fdf_get('SCF.MixAfterConvergence',compat_pre_v4_DM_H) 385 recompute_H_after_scf = fdf_get('SCF.Recompute-H-After-Scf',compat_pre_v4_DM_H) 386 387 if (ionode) then 388 if (compat_pre_v4_DM_H) then 389 write(6,2) ':!:Next two options activated by pre-4.0 compat. switch' 390 endif 391 write(6,1) 'redata: Mix DM or H after convergence',mix_after_convergence 392 write(6,1) 'redata: Recompute H after scf cycle', recompute_H_after_scf 393 endif 394 395 ! Pulay mixing, number of iterations for one Pulay mixing (maxsav) 396 maxsav = fdf_get('DM.NumberPulay', 0) 397 398 ! Broyden SCF mixing, number of iterations 399 broyden_maxit = fdf_get('DM.NumberBroyden',0) 400 ! FIRE SCF mixing, no parameters 401 fire_mix = fdf_get('DM.FIRE.Mixing',.false.) 402 403 if (cml_p) then 404 call cmlAddParameter( xf=mainXML, name='DM.NumberPulay', & 405 value=maxsav, dictRef='siesta:maxsav', & 406 units="cmlUnits:countable" ) 407 call cmlAddParameter( xf=mainXML, name='DM.NumberBroyden', & 408 value=broyden_maxit, dictRef='siesta:broyden_maxit', & 409 units="cmlUnits:countable" ) 410 endif 411 412 ! Mix density matrix on first SCF step (mix) 413 mix_scf_first = fdf_get('DM.MixSCF1', & 414 .not. compat_pre_v4_DM_H) 415 mix_scf_first = fdf_get('SCF.Mix.First', mix_scf_first) 416 mix_scf_first_force = fdf_get('SCF.Mix.First.Force', .false.) 417 if ( mix_scf_first_force ) then 418 ! Also set this, to note the user of mixing first SCF regardless 419 ! of flag. 420 mix_scf_first = .true. 421 end if 422 if (ionode) then 423 write(6,1) 'redata: Mix DM in first SCF step',mix_scf_first 424 endif 425 426 if (cml_p) then 427 call cmlAddParameter( xf=mainXML, name='DM.MixSCF1', & 428 value=mix_scf_first, dictRef='siesta:mix' ) 429 endif 430 431 ! Use disk or memory to store intermediate Pulay mixing vectors 432 ! (pulfile) 433 pulfile = fdf_get('DM.PulayOnFile',.false.) 434 if (ionode) then 435 if (pulfile) then 436 call die( 'redata: Cannot use DM.PulayOnFile=.true.'//& 437 'in this version' ) 438 endif 439 write(6,1) 'redata: Write Pulay info on disk',pulfile 440 endif 441 if (cml_p) then 442 call cmlAddParameter(xf=mainXML, name='DM.PulayOnFile', & 443 value=pulfile, dictRef='siesta:pulfile') 444 endif 445 446 ! Density Matrix Mixing (proportion of output DM in new input DM) 447 wmix = fdf_get('DM.MixingWeight',0.25_dp) 448!!$ if (ionode) then 449!!$ write(6,6) 'redata: New DM Mixing Weight',wmix 450!!$ endif 451 452 if (cml_p) then 453 call cmlAddParameter( xf=mainXML,name='DM.MixingWeight', & 454 value=wmix, dictRef='siesta:wmix', & 455 units="cmlUnits:dimensionless" ) 456 endif 457 458 ! Density Matrix occupancy tolerance 459 occtol = fdf_get('DM.OccupancyTolerance',1.e-12_dp) 460 if (ionode) then 461 write(6,8) 'redata: New DM Occupancy tolerance',occtol 462 endif 463 464 if (cml_p) then 465 call cmlAddParameter( xf=mainXML,name='DM.OccupancyTolerance', & 466 value=occtol, dictRef='siesta:occtol' , & 467 units="cmlUnits:dimensionless" ) 468 endif 469 470 ! Perform linear mixing each nkick SCF iterations (to kick system 471 ! when it is pinned in a poorly convergent SCF loop) 472 nkick = fdf_get('DM.NumberKick',0) 473 if (ionode) then 474 if (nkick .ge. 1) then 475 write(6,5) 'redata: Kick with linear mixing every',nkick,& 476 ' iterations' 477 else 478 write(6,2)'redata: No kicks to SCF' 479 endif 480 endif 481 482 if (cml_p) then 483 call cmlAddParameter( xf=mainXML,name='DM.NumberKick', & 484 value=nkick, dictRef='siesta:nkick', & 485 units="cmlUnits:countable" ) 486 endif 487 488 ! Density Matrix Mixing each nkick SCF iterations 489 wmixkick = fdf_get('DM.KickMixingWeight',0.5_dp) 490 if (ionode) then 491 write(6,6) 'redata: DM Mixing Weight for Kicks',wmixkick 492 endif 493 494 if (cml_p) then 495 call cmlAddParameter( xf=mainXML,name='DM.KickMixingWeight', & 496 value=wmixkick, dictRef='siesta:wmixkick',& 497 units="cmlUnits:dimensionless" ) 498 endif 499 500 501 !-----------------------------------! 502 ! ! 503 ! Reading of convergence criteria ! 504 ! ! 505 !-----------------------------------! 506 507 ! Harris energy convergence 508 ! If this is TRUE, then the following options are defaulted to FALSE. 509 converge_Eharr = fdf_get('DM.RequireHarrisConvergence', .false.) 510 converge_Eharr = fdf_get('SCF.Harris.Converge', converge_Eharr) 511 tolerance_Eharr = fdf_get('DM.HarrisTolerance', 1.e-4_dp*eV, 'Ry' ) 512 tolerance_Eharr = fdf_get('SCF.Harris.Tolerance', tolerance_Eharr, 'Ry' ) 513 if ( IONode ) then 514 write(6,1) 'redata: Require Harris convergence for SCF', & 515 converge_Eharr 516 write(6,7) 'redata: Harris energy tolerance for SCF', tolerance_Eharr/eV, ' eV' 517 end if 518 519 if (cml_p) then 520 call cmlAddParameter( xf=mainXML, name='SCF.Harris.Converge', & 521 value=converge_Eharr, & 522 dictRef='siesta:ReqHarrisConv' ) 523 call cmlAddParameter( xf=mainXML, name='SCF.Harris.Tolerance', units='siestaUnits:eV', & 524 value=tolerance_Eharr/eV, dictRef='siesta:Harris_tolerance') 525 end if 526 527 528 ! Density matrix convergence 529 converge_DM = fdf_get('SCF.DM.Converge', .not. converge_Eharr) 530 dDtol = fdf_get('DM.Tolerance',1.e-4_dp) 531 dDtol = fdf_get('SCF.DM.Tolerance',dDtol) 532 if ( IONode ) then 533 write(6,1) 'redata: Require DM convergence for SCF', converge_DM 534 write(6,11) 'redata: DM tolerance for SCF',dDtol 535 end if 536 if (cml_p) then 537 call cmlAddParameter( xf=mainXML, name='SCF.DM.Converge', & 538 value=converge_DM, & 539 dictRef='siesta:ReqDMConv' ) 540 call cmlAddParameter( xf=mainXML, name='SCF.DM.Tolerance', & 541 value=dDtol, dictRef='siesta:dDtol', & 542 units='siestaUnits:eAng_3' ) 543 end if 544 545 546 ! Energy-density matrix convergence 547 converge_EDM = fdf_get('SCF.EDM.Converge', .false.) 548 tolerance_EDM = fdf_get('SCF.EDM.Tolerance',1.e-3_dp*eV, 'Ry') 549 if ( IONode ) then 550 write(6,1) 'redata: Require EDM convergence for SCF', converge_EDM 551 write(6,7) 'redata: EDM tolerance for SCF',tolerance_EDM/eV, ' eV' 552 end if 553 if (cml_p) then 554 call cmlAddParameter( xf=mainXML, name='SCF.EDM.Converge', & 555 value=converge_EDM, & 556 dictRef='siesta:ReqEDMConv' ) 557 call cmlAddParameter( xf=mainXML, name='SCF.EDM.Tolerance', & 558 value=tolerance_EDM/eV, dictRef='siesta:EDM_tolerance', & 559 units='siestaUnits:eVeAng_3' ) 560 end if 561 562 563 ! Hamiltonian convergence 564 if ( compat_pre_v4_DM_H ) then 565 tBool = .false. 566 else 567 tBool = .not. converge_Eharr 568 end if 569 converge_H = fdf_get('SCF.H.Converge', tBool) 570 dHtol = fdf_get('SCF.H.Tolerance',1.e-3_dp*eV, 'Ry') 571 if ( IONode ) then 572 write(6,1) 'redata: Require H convergence for SCF', converge_H 573 write(6,7) 'redata: Hamiltonian tolerance for SCF', dHtol/eV, ' eV' 574 end if 575 576 if (cml_p) then 577 call cmlAddParameter( xf=mainXML, name='SCF.H.Converge', & 578 value=converge_H, & 579 dictRef='siesta:ReqHConv' ) 580 call cmlAddParameter( xf=mainXML, name='SCF.H.Tolerance', & 581 value=dHtol/eV, dictRef='siesta:dHtol', & 582 units='siestaUnits:eV' ) 583 end if 584 585 586 ! Free energy convergence 587 converge_FreeE = fdf_get('DM.RequireEnergyConvergence', .false.) 588 converge_FreeE = fdf_get('SCF.FreeE.Converge',converge_FreeE) 589 tolerance_FreeE = fdf_get('DM.EnergyTolerance', 1.e-4_dp*eV, 'Ry' ) 590 tolerance_FreeE = fdf_get('SCF.FreeE.Tolerance',tolerance_FreeE, 'Ry') 591 if ( IONode ) then 592 write(6,1) 'redata: Require (free) Energy convergence for SCF', converge_FreeE 593 write(6,7) 'redata: (free) Energy tolerance for SCF', tolerance_FreeE/eV, ' eV' 594 end if 595 596 if (cml_p) then 597 call cmlAddParameter( xf=mainXML, name='SCF.FreeE.Converge', & 598 value=converge_FreeE, & 599 dictRef='siesta:ReqEnergyConv' ) 600 call cmlAddParameter( xf=mainXML, name='SCF.FreeE.Tolerance', & 601 value=tolerance_FreeE/eV, dictRef='siesta:dEtol', & 602 units="siestaUnits:eV" ) 603 end if 604 605 ! Check that there indeed is at least one convergence criteria 606 tBool = .false. 607 tBool = tBool .or. converge_Eharr 608 tBool = tBool .or. converge_FreeE 609 tBool = tBool .or. converge_EDM 610 tBool = tBool .or. converge_DM 611 tBool = tBool .or. converge_H 612 if ( .not. tBool ) then 613 call die('There is no convergence criteria. Please at least have one.') 614 end if 615 616 !------------------------------------ 617 ! Done reading convergence criteria 618 !------------------------------------ 619 620 ! Monitor forces and stresses during SCF loop 621 monitor_forces_in_scf = fdf_get('MonitorForcesInSCF',.false.) 622 monitor_forces_in_scf = fdf_get('SCF.MonitorForces',monitor_forces_in_scf) 623 624 !-------------------------------------- 625 ! Initial spin density: Maximum polarization, Ferro (false), AF (true) 626 if ( nspin .eq. 2 ) then 627 init_anti_ferro = fdf_get('DM.InitSpin.AF',.false.) 628 if ( ionode ) then 629 write(6,1) 'redata: Antiferro initial spin density',init_anti_ferro 630 endif 631 if (cml_p) then 632 call cmlAddParameter( xf=mainXML, name='DM.InitSpinAF', & 633 value=init_anti_ferro, dictRef='siesta:inspn') 634 end if 635 end if 636 637 ! Use Saved Data 638 usesaveddata = fdf_get('UseSaveData',.false.) 639 if (ionode) then 640 write(6,1) 'redata: Using Saved Data (generic)', usesaveddata 641 endif 642 643 ! Use continuation DM files 644 usesavedm = fdf_get('DM.UseSaveDM',usesaveddata) 645 if (ionode) then 646 write(6,1) 'redata: Use continuation files for DM', usesavedm 647 endif 648 649 if (cml_p) then 650 call cmlAddParameter( xf=mainXML, name='DM.UseSaveDM', & 651 value=usesavedm, dictRef='siesta:usesavedm') 652 endif 653 654 ! Neglect Interactions between non-overlapping orbitals ... 655 negl = fdf_get('NeglNonOverlapInt',.false.) 656 if (ionode) then 657 write(6,1) 'redata: Neglect nonoverlap interactions',negl 658 endif 659 660 if (cml_p) then 661 call cmlAddParameter( xf=mainXML, name='NeglNonOverlapInt', & 662 value=negl, dictRef='siesta:negl' ) 663 endif 664 665 ! Method to Solve LDA Hamiltonian ... 666 method = fdf_get('SolutionMethod','diagon') 667 if (cml_p) then 668 call cmlAddParameter( xf=mainXML, name='SolutionMethod', & 669 value=method, dictRef='siesta:SCFmethod' ) 670 endif 671 672 if (leqi(method,'matrix')) then 673 isolve = MATRIX_WRITE 674 if (ionode) then 675 write(*,3) 'redata: Method of Calculation', 'Matrix write only' 676 endif 677 678 else if (leqi(method,'diagon')) then 679 isolve = SOLVE_DIAGON 680 if (ionode) then 681 write(*,3) 'redata: Method of Calculation', 'Diagonalization' 682 endif 683 684 else if (leqi(method,'ordern')) then 685 isolve = SOLVE_ORDERN 686 if (ionode) then 687 write(*,3) 'redata: Method of Calculation','Order-N' 688 endif 689 if (nspin .gt. 2) then 690 call die( 'redata: You chose the Order-N solution option '// & 691 'together with nspin>2. This is not allowed in '//& 692 'this version of siesta' ) 693 endif 694 695 else if (leqi(method,'omm')) then 696 isolve = SOLVE_MINIM 697 call_diagon_default=fdf_integer('OMM.Diagon',0) 698 call_diagon_first_step=fdf_integer('OMM.DiagonFirstStep',call_diagon_default) 699 minim_calc_eigenvalues=fdf_boolean('OMM.Eigenvalues',.false.) 700 if (ionode) then 701 write(*,3) 'redata: Method of Calculation', 'Orbital Minimization Method' 702 endif 703 else if (leqi(method,"pexsi")) then 704#ifdef SIESTA__PEXSI 705 isolve = SOLVE_PEXSI 706 if (ionode) then 707 call add_citation("10.1088/0953-8984/26/30/305503") 708 write(*,3) 'redata: Method of Calculation', 'PEXSI' 709 endif 710#else 711 call die("PEXSI solver is not compiled in. Use -DSIESTA__PEXSI") 712#endif 713 714 else if (leqi(method,'transi') .or. leqi(method,'transiesta') & 715 .or. leqi(method,'negf') ) then 716 isolve = SOLVE_TRANSI 717 if (ionode) then 718 call add_citation("10.1103/PhysRevB.65.165401") 719 call add_citation("10.1016/j.cpc.2016.09.022") 720 write(*,3) 'redata: Method of Calculation','Transiesta' 721 endif 722 if ( nspin > 2 ) then 723 call die('transiesta does not work for non-collinear or spin-orbit') 724 end if 725 else 726 call die( 'redata: The method of solution must be either '//& 727 'Transiesta, '//& 728#ifdef SIESTA__PEXSI 729 'PEXSI, '//& 730#endif 731 'OrderN, OMM or Diagon' ) 732 endif 733 734#ifdef DEBUG 735 call write_debug( ' Solution Method: ' // method ) 736#endif 737 738 ! Electronic temperature for Fermi Smearing ... 739 temp = fdf_get('ElectronicTemperature',1.9e-3_dp,'Ry') 740 if (ionode .and. isolve.eq.SOLVE_DIAGON) then 741 write(6,6) 'redata: Electronic Temperature',temp/Kelvin,' K' 742 endif 743 744 if (cml_p) then 745 call cmlAddParameter( xf=mainXML, name='ElectronicTemperature', & 746 value=temp, dictRef='siesta:etemp', & 747 units='siestaUnits:Ry') 748 endif 749 750 ! Fix the spin of the system to a given value ; and 751 ! value of the Spin of the system (only used if fixspin = TRUE) 752 fixspin = fdf_get('FixSpin',.false.) 753 fixspin = fdf_get('Spin.Fix',fixspin) 754 if (ionode) then 755 write(6,1) 'redata: Fix the spin of the system',fixspin 756 endif 757 758 if (fixspin) then 759 if (nspin .ne. 2) then 760 call die( 'redata: ERROR: You can only fix the spin of '//& 761 'the system for collinear spin polarized calculations.' ) 762 endif 763 total_spin = fdf_get('TotalSpin',0.0_dp) 764 total_spin = fdf_get('Spin.Total',total_spin) 765 if (ionode) then 766 write(6,9) 'redata: Total spin of the system (spin value)', total_spin 767 endif 768 else 769 total_spin = 0.0_dp 770 endif 771 772 if (cml_p) then 773 call cmlAddParameter( xf=mainXML, name='FixSpin', & 774 value=fixspin, dictref='siesta:fixspin' ) 775 call cmlAddParameter( xf=mainXML, name='TotalSpin', & 776 value=total_spin, dictref='siesta:totalspin',& 777 units='siestaUnits:eSpin' ) 778 endif 779 780 ! Order-N solution parameters ... 781 ! Maximum number of CG minimization iterations 782 ncgmax = fdf_get('ON.MaxNumIter',1000) 783 if (ncgmax<1) then 784 if (ionode) then 785 write(6,2) 'ON.MaxNumIter cannot be less than 1. Resetting to 1' 786 endif 787 ncgmax = 1 788 endif 789 790 ! Relative tolerance in total band structure energy 791 etol = fdf_get('ON.etol',1.e-8_dp) 792 793 ! Fermi level parameter 794 eta(1:2) = 0.0_dp 795 eta(1) = fdf_physical('ON.eta',eta(1),'Ry') 796 eta(2) = eta(1) 797 eta(1) = fdf_physical('ON.eta_alpha',eta(1),'Ry') 798 eta(2) = fdf_physical('ON.eta_beta',eta(2),'Ry') 799 800 ! Cutoff radius for Localized Wave Functions 801 rcoor = fdf_get('On.RcLWF',9.5_dp,'Bohr') 802 803 ! Use continumation LWF files 804 usesavelwf = fdf_get('ON.UseSaveLWF',usesaveddata) 805 806 ! Option on how to build LWF's (disk or functionals) 807 lwfopt = fdf_get('ON.functional','kim') 808 if (leqi(lwfopt,'files')) then 809 ioptlwf = 0 810 else if (leqi(lwfopt,'kim')) then 811 ioptlwf = 1 812 else if (leqi(lwfopt,'ordejon-mauri')) then 813 ioptlwf = 2 814 else 815 call die('redata: wrong ON.funcional option') 816 endif 817 818 ! Option to calculate the Chemical potential in O(N) 819 ! Option to use the Chemical Potential calculated instead 820 ! of the eta variable of the input 821 noeta = fdf_get('ON.ChemicalPotentialUse',.false.) 822 ! NOTE: This does not yet work in parallel 823 824 if (noeta) then 825 ! if so, we must (obviously) calculate the chemical potential 826 chebef=.true. 827 else 828 ! otherwise, we may still want to calculate it but not use it. 829 chebef = fdf_get('ON.ChemicalPotential',.false.) 830 endif 831 832#ifdef MPI 833 if (chebef) then 834 call die("ON.ChemicalPotential(Use) options do not work with MPI") 835 endif 836#endif 837 838 ! Cutoff radius to calculate the Chemical Potential by projection 839 rcoorcp = fdf_get( 'ON.ChemicalPotentialRc', 9.5_dp, 'Bohr' ) 840 841 ! Temperature of the Fermi distribution to calculate the 842 ! Chemical potential by projection 843 tcp = fdf_get( 'ON.ChemicalPotentialTemperature', 0.05_dp,'Ry' ) 844 beta = 1.0_dp/tcp 845 846 ! Order of the Chebishev expansion to calculate the Chemical potential 847 pmax = fdf_get('ON.ChemicalPotentialOrder',100) 848 849 850 if (isolve==SOLVE_ORDERN) then 851 if (ionode) then 852 write(6,4) 'redata: Maximum number of iterations',ncgmax 853 write(6,6) 'redata: Relative tolerance',etol 854 if (nspin.eq.2) then 855 write(6,6) 'redata: Eta (Fermi level) Alpha spin',eta(1)/eV,' eV' 856 write(6,6) 'redata: Eta (Fermi level) Beta spin',eta(2)/eV,' eV' 857 else 858 write(6,6) 'redata: Eta (Fermi level parameter)',eta(1)/eV,' eV' 859 endif 860 write(6,6) 'redata: Radius of LWFs',rcoor/Ang,' Ang' 861 write(6,1) 'redata: Use continuation files for LWF',usesavelwf 862 write(6,3) 'redata: Method to build LWFs',lwfopt 863 if (chebef) then 864 write(6,1) 'redata: Compute Chemical Potential',chebef 865 write(6,2) 'redata: Use the calculated Chemical ..' 866 write(6,1) 'redata: ..Potential instead of eta',noeta 867 write(6,6) 'redata: Radius to compute the Chem. Pot.',rcoorcp/Ang,' Ang' 868 write(6,2) 'redata: Temp. for Fermi distribution ..' 869 write(6,6) 'redata: .. to compute the Chem. Pot.',tcp/eV,' eV' 870 write(6,4) 'redata: Order of the Chebishev expansion',pmax 871 endif 872 endif 873 if (cml_p) then 874 call cmlAddParameter( xf = mainXML, & 875 name = 'ON.MaxNumIter',& 876 value = ncgmax, & 877 dictref = 'siesta:ncgmax', & 878 units = "cmlUnits:countable" ) 879 880 call cmlAddParameter( xf = mainXML, & 881 name = 'ON.etol', & 882 value = etol, & 883 dictref = 'siesta:etol', & 884 units = "siestaUnits:eV" ) 885 if (nspin==2) then 886 call cmlAddParameter( xf = mainXML, & 887 name = 'ON.eta_alpha', & 888 value = eta(1), & 889 dictref = 'siesta:eta1', & 890 units = 'siestaUnits:Ry' ) 891 892 call cmlAddParameter( xf = mainXML, & 893 name = 'ON.eta_beta', & 894 value = eta(2), & 895 dictref = 'siesta:eta2', & 896 units = 'siestaUnits:Ry' ) 897 else 898 call cmlAddParameter( xf = mainXML, & 899 name = 'ON.eta', & 900 value = eta(1), & 901 dictref = 'siesta:eta', & 902 units = 'siestaUnits:Ry') 903 endif 904 call cmlAddParameter( xf = mainXML, & 905 name = 'On.RcLWF', & 906 value = rcoor, & 907 dictref = 'siesta:rcoor', & 908 units = 'siestaUnits:Bohr') 909 910 call cmlAddParameter( xf=mainXML, & 911 name='On.UseSaveLWF', & 912 value=usesavelwf, & 913 dictref='siesta:usesavelwf' ) 914 915 call cmlAddParameter( xf = mainXML, & 916 name = 'ON.functional',& 917 value = lwfopt, & 918 dictref = 'siesta:lwfopt' ) 919 if (chebef) then 920 call cmlAddParameter( xf = mainXML, & 921 name = 'ON.ChemicalPotential', & 922 value = chebef, & 923 dictref = 'siesta:chebef') 924 925 call cmlAddParameter( xf = mainXML, & 926 name = 'ON.ChemicalPotentialUse', & 927 value = noeta, & 928 dictref = 'siesta:noeta') 929 930 call cmlAddParameter( xf = mainXML, & 931 name = 'ON.ChemicalPotentialRc', & 932 value = rcoorcp, & 933 dictref = 'siesta:rcoorcp', & 934 units = 'siestaUnits:Bohr') 935 936 call cmlAddParameter( xf = mainXML, & 937 name = 'ON.ChemicalPotentialTemperature', & 938 value = tcp, dictref='siesta:tcp', & 939 units = 'siestaUnits:Ry' ) 940 941 call cmlAddParameter( xf = mainXML, & 942 name = 'ON.ChemicalPotentialOrder', & 943 value = pmax, & 944 dictref = 'siesta:pmax', & 945 units = 'cmlUnits:dimensionless') 946 endif 947 endif 948 endif 949 950 ! Dynamics parameters ... 951 varcel = fdf_get('MD.VariableCell', .false. ) 952 953 ! NB reset below ... 954 ! Type of dynamics 955 956 compat_pre_v4_dynamics = fdf_get('compat-pre-v4-dynamics', .false. ) 957 if (compat_pre_v4_dynamics) then 958 dyntyp = fdf_get('MD.TypeOfRun','verlet') 959 else 960 dyntyp = fdf_get('MD.TypeOfRun','cg') 961 endif 962 963 if (leqi(dyntyp,'cg')) then 964 idyn = 0 965 usesavecg = fdf_get('MD.UseSaveCG', usesaveddata) 966 ! Support the old Broyden switch for now 967 broyden_optim = fdf_get('Optim.Broyden',.false.) 968 call deprecated('Optim.Broyden') 969 970 else if (leqi(dyntyp,'broyden')) then 971 idyn = 0 972 broyden_optim = .true. 973 else if (leqi(dyntyp,'fire')) then 974 idyn = 0 975 fire_optim = .true. 976 else if (leqi(dyntyp,'verlet')) then 977 idyn = 1 978 else if (leqi(dyntyp,'nose')) then 979 idyn = 2 980 else if (leqi(dyntyp,'parrinellorahman')) then 981 idyn = 3 982 else if (leqi(dyntyp,'noseparrinellorahman')) then 983 idyn = 4 984 else if (leqi(dyntyp,'anneal')) then 985 idyn = 5 986 else if (leqi(dyntyp,'fc')) then 987 idyn = 6 988 else if (leqi(dyntyp,'phonon')) then 989 call die('Dynamics type "PHONON" is no longer supported') 990 else if (leqi(dyntyp,'forces').or.leqi(dyntyp,'master')) then 991 idyn = 8 992#ifdef NCDF_4 993 else if (leqi(dyntyp,'explicit')) then 994 idyn = 9 995#endif 996#ifdef SIESTA__FLOOK 997 else if (leqi(dyntyp,'lua')) then 998 idyn = 10 999#endif 1000 else 1001 call die('Invalid Option selected - value of MD.TypeOfRun not recognised') 1002 endif 1003 1004 ! Maximum number of steps in MD/coordinate optimization 1005 nmove = fdf_get('MD.NumCGsteps',0) 1006 nmove = fdf_get('MD.Steps',nmove) 1007 1008 ! Maximum atomic displacement in one step 1009 dxmax = fdf_get('MD.MaxCGDispl',0.2_dp,'Bohr') 1010 dxmax = fdf_get('MD.MaxDispl',dxmax,'Bohr') 1011 1012 ! Tolerance in the maximum atomic force [0.04 eV/Ang] 1013 ftol = fdf_get('MD.MaxForceTol', 0.00155574_dp, 'Ry/Bohr') 1014 1015 ! Tolerance in the maximum residual stress (var cell) [1 GPa] 1016 strtol = fdf_get('MD.MaxStressTol', 6.79773e-5_dp, 'Ry/Bohr**3') 1017 strtol = abs(strtol) 1018 1019 GeometryMustConverge = fdf_get('GeometryMustConverge', .false.) 1020 1021 if (ionode) then 1022 select case (idyn) 1023 case(0) 1024 if (nmove > 0) then 1025 if (broyden_optim) then 1026 write(6,3) 'redata: Dynamics option','Broyden coord. optimization' 1027 elseif (fire_optim) then 1028 write(6,3) 'redata: Dynamics option', 'FIRE coord. optimization' 1029 else 1030 write(6,3) 'redata: Dynamics option','CG coord. optimization' 1031 endif 1032 write(6,1) 'redata: Variable cell', varcel 1033 if (.not. broyden_optim) then 1034 write(6,1) 'redata: Use continuation files for CG', usesavecg 1035 write(6,6) 'redata: Max atomic displ per move', dxmax/Ang, ' Ang' 1036 endif 1037 write(6,4) 'redata: Maximum number of optimization moves', nmove 1038 write(6,6) 'redata: Force tolerance', ftol/eV*Ang, ' eV/Ang' 1039 if (varcel) then 1040 write(6,6) 'redata: Stress tolerance', & 1041 strtol/6.79773e-5_dp, ' GPa' 1042 endif 1043 if (cml_p) then 1044 if (broyden_optim) then 1045 call cmlAddParameter( xf = mainXML, & 1046 name = 'MD.TypeOfRun', & 1047 value= 'Broyden' ) 1048 else if (fire_optim) then 1049 call cmlAddParameter( xf = mainXML, & 1050 name = 'MD.TypeOfRun', & 1051 value= 'FIRE' ) 1052 else 1053 call cmlAddParameter( xf =mainXML, & 1054 name ='MD.TypeOfRun', & 1055 value ='CG' ) 1056 1057 call cmlAddParameter( xf = mainXML, & 1058 name = 'MD.UseSaveCG',& 1059 value = usesavecg ) 1060 endif 1061 1062 call cmlAddParameter( xf = mainXML, & 1063 name = 'MD.NumCGSteps', & 1064 value = nmove, & 1065 units = "cmlUnits:countable" ) 1066 call cmlAddParameter( xf = mainXML, & 1067 name = 'MD.Steps', & 1068 value = nmove, & 1069 units = "cmlUnits:countable" ) 1070 1071 call cmlAddParameter( xf = mainXML, & 1072 name = 'MD.MaxCGDispl', & 1073 value = dxmax, & 1074 units = 'siestaUnits:Bohr' ) 1075 call cmlAddParameter( xf = mainXML, & 1076 name = 'MD.MaxDispl', & 1077 value = dxmax, & 1078 units = 'siestaUnits:Bohr' ) 1079 1080 call cmlAddParameter( xf=mainXML, & 1081 name='MD.MaxForceTol', & 1082 value=ftol, & 1083 units='siestaUnits:Ry_Bohr') 1084 if (varcel) then 1085 call cmlAddParameter( xf=mainXML, & 1086 name='MD.MaxStressTol', & 1087 value=strtol, & 1088 units='siestaUnits:Ry_Bohr__3' ) 1089 endif 1090 endif 1091 else 1092 write(6,3) 'redata: Dynamics option','Single-point calculation' 1093 if (cml_p) then 1094 call cmlAddParameter( xf = mainXML, & 1095 name = 'MD.TypeOfRun', & 1096 value= 'Single-Point' ) 1097 endif 1098 endif 1099 case(1) 1100 write(6,3) 'redata: Dynamics option', 'Verlet MD run' 1101 if (cml_p) then 1102 call cmlAddParameter( xf = mainXML, & 1103 name = 'MD.TypeOfRun', & 1104 value = 'Verlet' ) 1105 endif 1106 1107 case(2) 1108 write(6,3) 'redata: Dynamics option', 'Nose thermostat MD run' 1109 if (cml_p) then 1110 call cmlAddParameter( xf=mainXML, name='MD.TypeOfRun', value='Nose') 1111 endif 1112 1113 case(3) 1114 write(6,3) 'redata: Dynamics option', 'Parrinello-Rahman MD run' 1115 if (cml_p) then 1116 call cmlAddParameter( xf= mainXML, & 1117 name= 'MD.TypeOfRun', & 1118 value= 'Parrinello-Rahman' ) 1119 endif 1120 1121 case(4) 1122 write(6,3) 'redata: Dynamics option', 'Nose-Parrinello-Rahman MD run' 1123 if (cml_p) then 1124 call cmlAddParameter( xf = mainXML, & 1125 name = 'MD.TypeOfRun', & 1126 value = 'Nose-Parrinello-Rahman' ) 1127 endif 1128 1129 case(5) 1130 write(6,3) 'redata: Dynamics option', 'Annealing MD run' 1131 if (cml_p) then 1132 call cmlAddParameter( xf = mainXML, & 1133 name = 'MD.TypeOfRun', & 1134 value = 'Annealing' ) 1135 endif 1136 1137 case(6) 1138 write(6,3) 'redata: Dynamics option', 'Force Constants Matrix Calculation' 1139 if (cml_p) then 1140 call cmlAddParameter( xf = mainXML, & 1141 name = 'MD.TypeOfRun', & 1142 value = 'Force Constants' ) 1143 endif 1144 1145 case(7) 1146 ! deprecated 1147 1148 case(8) 1149 write(6,3) 'redata: Dynamics option','Force evaluation' 1150 if (cml_p) then 1151 call cmlAddParameter( xf = mainXML, & 1152 name = 'MD.TypeOfRun', & 1153 value = 'Force Evaluation' ) 1154 endif 1155#ifdef NCDF_4 1156 case(9) 1157 write(6,3) 'redata: Dynamics option','Explicit' 1158 if (cml_p) then 1159 call cmlAddParameter( xf = mainXML, & 1160 name = 'MD.TypeOfRun', & 1161 value = 'Explicit' ) 1162 endif 1163#endif 1164#ifdef SIESTA__FLOOK 1165 case(10) 1166 write(6,3) 'redata: Dynamics option','LUA' 1167 if (cml_p) then 1168 call cmlAddParameter( xf = mainXML, & 1169 name = 'MD.TypeOfRun', & 1170 value = 'LUA' ) 1171 endif 1172#endif 1173 end select 1174 endif 1175 1176 ! Initial and final time steps for MD 1177 istart = fdf_get('MD.InitialTimeStep',1) 1178 if ( fdf_defined('MD.Steps') ) then 1179 ifinal = fdf_get('MD.FinalTimeStep',max(1,nmove - istart + 1)) 1180 else 1181 ifinal = fdf_get('MD.FinalTimeStep',1) 1182 end if 1183 1184 ! Length of time step for MD 1185 dt = fdf_get('MD.LengthTimeStep',1._dp,'fs') 1186 1187 ! Quench Option 1188 qnch = fdf_get('MD.Quench',.false.) 1189 qnch2 = fdf_get('MD.FireQuench',.false.) 1190 if ((qnch .or. qnch2) .and. (idyn==2 .or. idyn==4)) then 1191 call die( 'redata: ERROR: You cannot quench and '//& 1192 'use a Nose thermostat simultaneously') 1193 endif 1194 1195 iquench = 0 1196 if (qnch) then 1197 iquench = 1 1198 endif 1199 if (qnch2) then 1200 iquench = 2 1201 endif 1202 1203 ! Initial Temperature of MD simulation 1204 ! (draws random velocities from the Maxwell-Boltzmann distribition 1205 ! at the given temperature) 1206 tempinit = fdf_get('MD.InitialTemperature',0.0_dp,'K') 1207 1208 if (idyn .ge. 1 .and. idyn .le. 5) then 1209 if (ionode) then 1210 write(6,4) 'redata: Initial MD time step',istart 1211 write(6,4) 'redata: Final MD time step',ifinal 1212 write(6,6) 'redata: Length of MD time step',dt,' fs' 1213 write(6,6) 'redata: Initial Temperature of MD run',tempinit,' K' 1214 if (idyn .ne. 5) then 1215 if (qnch2) then 1216 write(6,1) 'redata: Perform a MD Fire quench',qnch2 1217 else 1218 write(6,1) 'redata: Perform a MD quench',qnch 1219 endif 1220 endif 1221 endif 1222 1223 if (cml_p) then 1224 call cmlAddParameter( xf = mainXML, & 1225 name = 'MD.InitialTimeStep',& 1226 value = istart, & 1227 units = 'cmlUnits:countable' ) 1228 1229 call cmlAddParameter( xf = mainXML, & 1230 name = 'MD.FinalTimeStep', & 1231 value = ifinal, & 1232 units = 'cmlUnits:countable' ) 1233 1234 call cmlAddParameter( xf=mainXML, & 1235 name='MD.LengthTimeStep',& 1236 value=dt, & 1237 units='siestaUnits:fs' ) 1238 1239 call cmlAddParameter( xf=mainXML, & 1240 name='MD.InitialTemperature',& 1241 value=tempinit, & 1242 units='siestaUnits:K' ) 1243 if (idyn/=5) then 1244 if(qnch2) then 1245 call cmlAddParameter( xf = mainXML, & 1246 name = 'MD.FireQuench', & 1247 value = qnch2 ) 1248 else 1249 call cmlAddParameter( xf = mainXML, & 1250 name = 'MD.Quench', & 1251 value = qnch ) 1252 endif 1253 endif 1254 endif 1255 endif 1256 1257 ! Target Temperature and Pressure 1258 tt = fdf_get('MD.TargetTemperature',0.0_dp,'K') 1259 tp = fdf_get('MD.TargetPressure',0.0_dp,'Ry/Bohr**3') 1260 ! 1261 ! Used for now for the call of the PR md routine if quenching 1262 if (idyn == 3 .AND. iquench > 0) call set_target_stress() 1263 1264 1265 ! Mass of Nose variable 1266 mn = fdf_get('MD.NoseMass',100._dp,'Ry*fs**2') 1267 1268 ! Mass of Parrinello-Rahman variables 1269 mpr = fdf_get('MD.ParrinelloRahmanMass',100._dp,'Ry*fs**2') 1270 1271 if (idyn==2 .or. idyn==4) then 1272 if (ionode) then 1273 write(6,6) 'redata: Nose mass',mn/eV,' eV/fs**2' 1274 endif 1275 if (cml_p) then 1276 call cmlAddParameter( xf = mainXML, & 1277 name = 'MD.NoseMass', & 1278 value = mn, & 1279 units = 'siestaUnits:Ry_fs__2') 1280 endif 1281 endif 1282 1283 if (idyn==3 .or. idyn==4) then 1284 if (ionode) then 1285 write(6,6) 'redata: Parrinello-Rahman mass',mpr/eV,' eV/fs**2' 1286 endif 1287 if (cml_p) then 1288 call cmlAddParameter( xf = mainXML, & 1289 name = 'MD.ParrinelloRahmanMass', & 1290 value = mn, & 1291 units = 'siestaUnits:Ry_fs__2' ) 1292 endif 1293 endif 1294 1295 ! Annealing option 1296 ianneal = 0 1297 annop = fdf_get( 'MD.AnnealOption','TemperatureAndPressure' ) 1298 1299 if (idyn .eq. 5) then 1300 if (leqi(annop,'Temperature')) then 1301 ianneal = 1 1302 else if (leqi(annop,'Pressure')) then 1303 ianneal = 2 1304 else if (leqi(annop,'TemperatureAndPressure')) then 1305 ianneal = 3 1306 else 1307 call die( 'redata: ERROR: With annealing MD, you must '//& 1308 'choose an appropriate value for MD.AnnealOption' ) 1309 endif 1310 1311 if (ionode) then 1312 select case (ianneal) 1313 case(1) 1314 write(6,3) 'redata: Annealing Option', 'Temperature' 1315 if (cml_p) then 1316 call cmlAddParameter( xf = mainXML, & 1317 name = 'MD.AnnealOption', & 1318 value = 'Temperature' ) 1319 endif 1320 1321 case(2) 1322 write(6,3) 'redata: Annealing Option', 'Pressure' 1323 if (cml_p) then 1324 call cmlAddParameter( xf = mainXML, & 1325 name = 'MD.AnnealOption', & 1326 value = 'Pressure') 1327 endif 1328 1329 case(3) 1330 write(6,3) 'redata: Annealing Option', 'Temperature and Pressure' 1331 if (cml_p) then 1332 call cmlAddParameter( xf = mainXML, & 1333 name = 'MD.AnnealOption', & 1334 value = 'TemperatureAndPressure') 1335 endif 1336 end select 1337 endif 1338 endif 1339 1340 1341 if (idyn==2 .or. idyn==4 .or. (idyn==5 .and. (ianneal ==1 .or. ianneal==3))) then 1342 if (ionode) then 1343 write(6,6) 'redata: Target Temperature',tt,' Kelvin' 1344 endif 1345 1346 if (cml_p) then 1347 call cmlAddParameter( xf = mainXML, & 1348 name = 'MD.TargetTemperature', & 1349 value = tt, & 1350 units = 'siestaUnits:K' ) 1351 endif 1352 endif 1353 1354 if (idyn==3 .or. idyn==4 .or. (idyn==5 .and. (ianneal==2 .or. ianneal==3))) then 1355 if (ionode) then 1356 write(6,6) 'redata: Target Pressure', tp/eV*Ang**3, ' eV/Ang**3' 1357 endif 1358 1359 if (cml_p) then 1360 call cmlAddParameter( xf=mainXML, & 1361 name= 'MD.TargetPressure', & 1362 value= tp, & 1363 units= 'siestaUnits:Ry_Bohr__3' ) 1364 endif 1365 endif 1366 1367 ! Relaxation Time for Annealing 1368 taurelax = fdf_get( 'MD.TauRelax', 100._dp,'fs' ) 1369 if (idyn==5) then 1370 if (ionode) then 1371 write(6,6) 'redata: Annealing Relaxation Time', taurelax,' fs' 1372 endif 1373 if (cml_p) then 1374 call cmlAddParameter( xf = mainXML, & 1375 name = 'MD.TauRelax', & 1376 value = taurelax, & 1377 units = 'siestaUnits:fs' ) 1378 endif 1379 endif 1380 1381 ! Estimated Bulk modulus (for Pressure annealing) [100 GPa] 1382 bulkm = fdf_get( 'MD.BulkModulus',100*6.79773e-5_dp,'Ry/Bohr**3' ) 1383 if (ionode) then 1384 if (idyn==5 .and. (ianneal==2 .or. ianneal==3)) then 1385 write(6,6) 'redata: Approx. Bulk Modulus ', bulkm/eV*Ang**3, ' eV/Ang**3' 1386 endif 1387 endif 1388 1389 if (cml_p) then 1390 call cmlAddParameter( xf = mainXML, & 1391 name = 'MD.BulkModulus', & 1392 value = bulkm, & 1393 units = 'siestaUnits:Ry_Bohr__3') 1394 endif 1395 1396 ! Atomic displacement for force constant calculation 1397 dx = fdf_get('MD.FCDispl',0.04_dp,'Bohr') 1398 1399 ! First and last atoms to displace for calculation of force constants 1400 ia1 = fdf_get('MD.FCfirst',1) 1401 ia2 = fdf_get('MD.FClast',na) 1402 1403 ! Check that last atom doesn't exceed total number 1404 if (idyn.eq.6.and.ia2.gt.na) then 1405 call die( 'redata: ERROR:'//& 1406 'Last atom index for FC calculation is > number of atoms.') 1407 endif 1408 1409 if (idyn==6) then 1410 if (ionode) then 1411 write(6,6) 'redata: Atomic displ for force constants',dx/Ang,' Ang' 1412 write(6,4) 'redata: First atom to move',ia1 1413 write(6,4) 'redata: Last atom to move',ia2 1414 endif 1415 if (cml_p) then 1416 call cmlAddParameter( xf = mainXML, & 1417 name = 'MD.FCDispl', & 1418 value = dx, & 1419 units = 'siestaUnits:Bohr' ) 1420 1421 call cmlAddParameter( xf= mainXML, & 1422 name= 'MD.FCFirst', & 1423 value= ia1, & 1424 units= 'cmlUnits:countable' ) 1425 1426 call cmlAddParameter( xf= mainXML, & 1427 name= 'MD.FCLast', & 1428 value= ia2, & 1429 units= 'cmlUnits:countable' ) 1430 1431 endif 1432 endif 1433 1434 ! Variable cell shape? Depending on input and type of dynamics 1435 varcel = varcel .or. (idyn==3) .or. (idyn==4) & 1436 .or. (idyn==5 .and. ianneal==1) & 1437 .and. (idyn/=1) .and. (idyn/=2) & 1438 .and. (idyn/=6) .and. (idyn/=7) & 1439 .and. (idyn/=9) & 1440 .and. (.not. (idyn==5 .and. ianneal/=1) ) 1441 1442 want_spatial_decomposition = fdf_get('UseSpatialDecomposition', .false.) 1443 want_domain_decomposition = fdf_get('UseDomainDecomposition', .false.) 1444#if defined(ON_DOMAIN_DECOMP) || defined(SIESTA__METIS) 1445#else 1446#ifdef MPI 1447 if (want_domain_decomposition) then 1448 write(*,*)'Need to compile SIESTA__METIS or ON_DOMAIN_DECOMP' 1449 call die("Need to compile with ON_DOMAIN_DECOMP support") 1450 endif 1451#endif 1452#endif 1453 1454 ! Read in mixing parameters (SCF) 1455 call mixers_scf_init( nspin ) 1456 call mixers_scf_print( nspin ) 1457 1458 ! We read in relevant data for ChargeGeometries block 1459 call read_charge_add( min(2, nspin) , charnet ) 1460 1461 ! We read in the relevant data for HartreeGeometries block 1462 call read_hartree_add( ) 1463 1464 ! Harris Forces?. Then DM.UseSaveDM should be false (use always 1465 ! Harris density in the first SCF step of each MD step), and 1466 ! MaxSCFIter should be 2, in the second one the Harris 1467 ! forces are computed. Also, should not exit if SCF did 1468 ! not converge. 1469 1470 harrisfun = fdf_get('Harris_functional',.false.) 1471 1472 ! First read in, then later correct depending on 1473 ! the other usages 1474 use_aux_cell = fdf_get('ForceAuxCell', .false.) 1475 1476 if (harrisfun) then 1477 mixH = .false. 1478 mix_charge = .false. 1479 usesavedm = .false. 1480 nscf = 1 ! Note change from tradition, since siesta_forces 1481 ! now explicitly separates the "compute_forces" 1482 ! phase from the rest of the scf cycle. 1483 mix_scf_first = .false. 1484 SCFMustConverge = .false. 1485 if ( nspin > 2 ) then 1486 call die("Harris functional does not work for & 1487 &non-collinear spin polarization.") 1488 end if 1489 endif 1490 1491 ! Warn the user about deprecated symbols... 1492 call deprecated('Optim.Broyden.History.Steps') 1493 call deprecated('Optim.Broyden.Cycle.On.Maxit') 1494 call deprecated('Optim.Broyden.Variable.Weight') 1495 call deprecated('Optim.Broyden.Debug') 1496 call deprecated('Optim.Broyden.Initial.Inverse.Jacobian') 1497 1498 1499 ! Find some switches 1500 writek = fdf_get( 'WriteKpoints', outlng ) 1501 writef = fdf_get( 'WriteForces', outlng ) 1502 1503 writedm = fdf_get( 'WriteDM', .true. ) 1504 write_dm_at_end_of_cycle = fdf_get( 'WriteDM.End.Of.Cycle', writedm ) 1505 writeH = fdf_get( 'WriteH', .false. ) 1506 write_H_at_end_of_cycle = fdf_get( 'WriteH.End.Of.Cycle', writeH ) 1507 1508 writedm_cdf = fdf_get('WriteDM.NetCDF', .false. ) 1509 writedm_cdf_history = fdf_get('WriteDM.History.NetCDF', .false. ) 1510 writedmhs_cdf = fdf_get('WriteDMHS.NetCDF', .false. ) 1511 writedmhs_cdf_history = fdf_get('WriteDMHS.History.NetCDF', .false.) 1512 read_charge_cdf = fdf_get('SCF.Read.Charge.NetCDF' , .false. ) 1513 read_deformation_charge_cdf = & 1514 fdf_get('SCF.Read.Deformation.Charge.NetCDF', .false. ) 1515 1516 ! Write the history 1517 write_tshs_history = fdf_get('Write.TSHS.History', .false.) 1518 if ( IONode.and.write_tshs_history ) & 1519 write(*,2) 'redata: Saves TSHS files in MD simulation' 1520 !write_hs_history = fdf_get('Write.HS.History', .false.) 1521 1522 if (read_charge_cdf .or. read_deformation_charge_cdf) then 1523 mix_scf_first = .false. 1524 endif 1525 1526 save_initial_charge_density = fdf_get( & 1527 'SaveInitialChargeDensity' , .false.) 1528 1529 analyze_charge_density_only = fdf_get( & 1530 'AnalyzeChargeDensityOnly' , .false.) 1531 1532 tBool = fdf_get( 'UseNewDiagk', .false. ) 1533 if ( tBool ) then 1534 ctmp = fdf_get('Diag.WFS.Cache', 'CDF') 1535 else 1536 ctmp = fdf_get('Diag.WFS.Cache', 'none') 1537 end if 1538 if ( leqi(ctmp, 'none') ) then 1539 diag_wfs_cache = 0 1540 else if ( leqi(ctmp, 'cdf') ) then 1541 diag_wfs_cache = 1 1542 else 1543 call die('redata: ERROR: Diag.WFS.Cache must be one of none|cdf') 1544 end if 1545 1546 writb = fdf_get( 'WriteBands', outlng ) 1547 writbk = fdf_get( 'WriteKbands', outlng ) 1548 writeig = fdf_get('WriteEigenvalues', outlng ) 1549 writec = fdf_get( 'WriteCoorStep', outlng ) 1550 writmd = fdf_get( 'WriteMDhistory', .false. ) 1551 writpx = fdf_get( 'WriteMDXmol', .not. writec ) 1552 save_ORB_INDX = fdf_get( 'WriteOrbitalIndex', .true. ) 1553 ! Do options of graphviz 1554 ctmp = fdf_get( 'Write.Graphviz', 'none' ) 1555 write_GRAPHVIZ = 0 1556 if ( leqi(ctmp, 'orb') .or. & 1557 leqi(ctmp, 'orbital') ) then 1558 write_GRAPHVIZ = 1 1559 else if ( leqi(ctmp, 'atom') .or. & 1560 leqi(ctmp, 'atomic') ) then 1561 write_GRAPHVIZ = 2 1562 else if ( leqi(ctmp, 'atom+orb') .or. & 1563 leqi(ctmp, 'atom+orbital') .or. & 1564 leqi(ctmp, 'orbital+atom') .or. & 1565 leqi(ctmp, 'orb+atom') .or. leqi(ctmp, 'both') ) then 1566 write_GRAPHVIZ = 3 1567 else if ( leqi(ctmp,'none') ) then 1568 ! do nothing, correct 1569 else if ( IONode ) then 1570 write(*,2) 'Write.Graphviz input could not be understood. & 1571 &Use [orbital|atom|atom+orbital] in the option.' 1572 end if 1573 if ( IONode ) then 1574 select case ( write_GRAPHVIZ ) 1575 case ( 1 ) 1576 write(*,2) 'redata: Save orbital connectivity graph in GRAPHVIZ format' 1577 case ( 2 ) 1578 write(*,2) 'redata: Save atomic connectivity graph in GRAPHVIZ format' 1579 case ( 3 ) 1580 write(*,2) 'redata: Save atom+orbital connectivity graphs in GRAPHVIZ format' 1581 end select 1582 end if 1583 1584 1585 writec = fdf_get( 'WriteCoorStep', outlng ) 1586 savehs = fdf_get( 'SaveHS', .false. ) 1587 initdmaux = fdf_get( 'ReInitialiseDM', .TRUE. ) 1588 allow_dm_reuse = fdf_get( 'DM.AllowReuse', .TRUE. ) 1589 allow_dm_extrapolation = fdf_get( 'DM.AllowExtrapolation', .TRUE. ) 1590 DM_history_depth = fdf_get( 'DM.History.Depth', 1) 1591 dm_normalization_tol = fdf_get( 'DM.NormalizationTolerance',1.0d-5) 1592 normalize_dm_during_scf= fdf_get( 'DM.NormalizeDuringSCF',.true.) 1593 muldeb = fdf_get( 'MullikenInSCF' , .false.) 1594 spndeb = fdf_get( 'SpinInSCF' , (nspin>1) ) 1595 1596 ! If no mulliken is requested, set it to false 1597 if ( mullipop == 0 ) muldeb = .false. 1598 rijmin = fdf_get( 'WarningMinimumAtomicDistance', & 1599 1.0_dp, 'Bohr' ) 1600 bornz = fdf_get( 'BornCharge' , .false. ) 1601 if (idyn.ne.6) then 1602 bornz = .false. 1603 endif 1604 change_kgrid_in_md = fdf_get('ChangeKgridInMD', .false.) 1605 RelaxCellOnly = fdf_get('MD.RelaxCellOnly', .false.) 1606 RemoveIntraMolecularPressure = fdf_get( & 1607 'MD.RemoveIntraMolecularPressure', .false.) 1608 ! 1609 ! COOP-related flags 1610 ! 1611 write_coop = fdf_get('COOP.Write', .false.) 1612 ! 1613 saverho = fdf_get( 'SaveRho', dumpcharge) 1614 savedrho = fdf_get( 'SaveDeltaRho', .false. ) 1615 saverhoxc= fdf_get( 'SaveRhoXC', .false.) 1616 savevh = fdf_get( 'SaveElectrostaticPotential', .false. ) 1617 savevna = fdf_get( 'SaveNeutralAtomPotential', .false. ) 1618 savevt = fdf_get( 'SaveTotalPotential', .false. ) 1619 savepsch = fdf_get( 'SaveIonicCharge', .false. ) 1620 savebader= fdf_get( 'SaveBaderCharge', .false.) 1621 savetoch = fdf_get( 'SaveTotalCharge', savebader ) 1622 1623 ! 1624 ! Siesta2Wannier90 -related flags 1625 ! 1626 w90_write_mmn = fdf_get( 'Siesta2Wannier90.WriteMmn', .false. ) 1627 w90_write_unk = fdf_get( 'Siesta2Wannier90.WriteUnk', .false. ) 1628 w90_write_amn = fdf_get( 'Siesta2Wannier90.WriteAmn', .false. ) 1629 w90_write_eig = fdf_get( 'Siesta2Wannier90.WriteEig', .false. ) 1630 1631 w90_processing = ( w90_write_mmn .or. w90_write_unk .or. & 1632 w90_write_amn .or. w90_write_eig ) 1633 1634 hasnobup = fdf_defined( 'Siesta2Wannier90.NumberOfBandsUp' ) 1635 hasnobdown = fdf_defined( 'Siesta2Wannier90.NumberOfBandsDown' ) 1636 hasnob = fdf_defined( 'Siesta2Wannier90.NumberOfBands' ) 1637 1638 nobup = fdf_get( 'Siesta2Wannier90.NumberOfBandsUp', 0) 1639 nobdown = fdf_get( 'Siesta2Wannier90.NumberOfBandsDown', 0) 1640 nob = fdf_get( 'Siesta2Wannier90.NumberOfBands', 0) 1641 1642#ifdef NCDF_4 1643 write_cdf = fdf_get('CDF.Save', .false.) 1644 ! No compression is by far the fastest 1645 cdf_comp_lvl = fdf_get('CDF.Compress', 0) 1646 if ( Nodes > 1 ) then 1647 cdf_w_parallel = fdf_get('CDF.MPI', .false.) 1648 else 1649 cdf_w_parallel = .false. 1650 end if 1651 1652 ! Only allow writing the CDF file for FC and non-MD calculations 1653 ! The MD file should be something different that only contains 1654 ! these things. 1655 if ( write_cdf ) then 1656 if ( idyn == 0 .and. nmove == 0 ) then 1657 ! non-MD calculation 1658 else if ( idyn == 6 ) then 1659 ! FC calculation, the FC calculation is fine 1660 ! Here we disable saving any real-space grid data 1661 save_initial_charge_density = .false. 1662 saverho = .false. 1663 savedrho = .false. 1664 saverhoxc = .false. 1665 savevh = .false. 1666 savevna = .false. 1667 savevt = .false. 1668 savepsch = .false. 1669 savebader = .false. 1670 savetoch = .false. 1671 else 1672 write_cdf = .false. 1673 end if 1674 end if 1675# ifndef NCDF_PARALLEL 1676 ! If not compiled with NCDF_PARALLEL, we do not 1677 ! allow parallel writes.....!!!! 1678 cdf_w_parallel = .false. 1679# endif 1680 if ( cdf_w_parallel ) then 1681 ! Doing parallel writes does not allow 1682 ! compression (the offsets cannot be calculated) 1683 cdf_comp_lvl = 0 1684 end if 1685 cdf_r_parallel = fdf_get('CDF.Read.Parallel', cdf_w_parallel ) 1686 1687 if ( IONode ) then 1688 ! Write out 1689 write(*,1) 'redata: Save all siesta data in one NC',write_cdf 1690 if ( write_cdf ) then 1691 if ( grid_p == dp ) then 1692 ctmp = fdf_get('CDF.Grid.Precision','double') 1693 if ( leqi(ctmp,'single') .or. leqi(ctmp,'float') ) then 1694 write(*,2) 'redata: Grids in NC reduced to single precision' 1695 end if 1696 end if 1697 write(*,4) 'redata: NC compression level',cdf_comp_lvl 1698 if ( cdf_r_parallel ) then 1699 write(*,2) 'redata: Reads NC in parallel' 1700 end if 1701 if ( cdf_w_parallel ) then 1702 write(*,2) 'redata: Writes NC in parallel (possibly not working)' 1703 end if 1704 end if 1705 end if 1706#endif 1707 1708 if (ionode) then 1709 write(6,'(2a)') 'redata: ', repeat('*', 71) 1710 endif 1711 1712 if (cml_p) then 1713 call cmlEndParameterList(mainXML) 1714 endif 1715 1716 ! Print blocks 1717 call mixers_scf_print_block( ) 1718 1719 RETURN 1720 !-------------------------------------------------------------------- END 17211 format(a,t53,'= ',2x,l1) 17222 format(a) 17233 format(a,t53,'= ',a) 17244 format(a,t53,'= ',i8) 17255 format(a,t53,'= ',i8,a) 17266 format(a,t53,'= ',f10.4,a) 17277 format(a,t53,'= ',f12.6,a) 17288 format(a,t53,'= ',f14.12) 17299 format(a,t53,'= ',f10.4) 173010 format(t55,a) 173111 format(a,t53,'= ',f12.6) 1732 1733CONTAINS 1734 subroutine deprecated( str ) 1735 1736 implicit none 1737 1738 character(len=*), intent(in) :: str 1739 1740 if (ionode) then 1741 if (fdf_defined(trim(str))) then 1742 write(6,'(a)') '**Warning: FDF symbol '//trim(str)//& 1743 ' is deprecated. See the manual.' 1744 endif 1745 endif 1746 1747 end subroutine deprecated 1748 1749end subroutine read_options 1750 1751 1752