1* 2* $Id$ 3* 4 5* *********************************** 6* * * 7* * control_read * 8* * * 9* *********************************** 10 logical function control_read(code_in,rtdb) 11 implicit none 12 integer code_in 13 integer rtdb 14 15#include "inp.fh" 16#include "bafdecls.fh" 17#include "btdb.fh" 18#include "control.fh" 19#include "nwpwxc.fh" 20#include "util.fh" 21 22 integer MASTER,taskid 23 parameter (MASTER=0) 24 25 logical value,np_default,value2,value3 26 integer ispin0,ne(2),nbrill,np1,np2 27 integer nmult(5) 28 real*8 error,thresh 29 30* **** control_nose common block **** 31 logical nose 32 real*8 Pe,Te,Pr,Tr 33 common / control_nblock / Pe,Te,Pr,Tr,nose 34 35* **** control_rtdb common block **** 36 integer trtdb 37 common / control_rtdb1 / trtdb 38 39* **** control_print common block **** 40 integer print_level 41 common / control_print1 / print_level 42 43 !character*30 cell_name 44 character*50 rtdb_unita,rtdb_unitaf,rtdb_ngrid,rtdb_boundry 45 character*50 exchange_correlation,rtdb_ngrid_small 46 integer i,l,np 47! integer ind ! unused 48 49 integer control_num_kvectors 50 external control_num_kvectors 51 logical control_allow_translation 52 external control_allow_translation 53 integer Parallel_threadid 54 external Parallel_threadid 55 56 57 value = .true. 58 59 call Parallel_taskid(taskid) 60 call nwpw_timing_start(50) 61 value2 = btdb_parallel(.true.) 62 code = code_in 63 trtdb = rtdb 64 65 call util_print_get_level(print_level) 66 !write(*,*) "Print level is ",print_level 67 68 69 DO 10, i = 1,100 70* **** get parallel mappings **** 71 if (.not.btdb_get(rtdb,'nwpw:mapping',mt_int,1,mapping)) then 72 mapping = 2 73 end if 7410 continue 75 76* **** set mapping1d **** 77 if (.not.btdb_get(rtdb,'nwpw:mapping1d',mt_int,1,mapping1d)) then 78 mapping1d = 1 79 end if 80 81* **** set np_dimensions **** 82 call Parallel_np(np) 83 84**** if NOMPI then no processor groups or communicators **** 85#ifdef NOMPI 86 np_dimensions(1) = np 87 np_dimensions(2) = 1 88 np_dimensions(3) = 1 89 np_default = .true. 90#else 91 if(.not.btdb_get(rtdb,'nwpw:np_dimensions', 92 > mt_int,3,np_dimensions)) then 93 94 np_dimensions(1) = np 95 np_dimensions(2) = 1 96 np_dimensions(3) = 1 97 np_default = .true. 98 else 99 np_default = .false. 100 if (.not.((code.eq.5) .or. 101 > (code.eq.10).or. 102 > (code.eq.13).or. 103 > (code.eq.14))) np_dimensions(3) = 1 104 if (np_dimensions(3).lt.1) np_dimensions(3) = 1 105 if (np_dimensions(2).lt.1) np_dimensions(2) = 1 106 if (np_dimensions(1).lt.1) np_dimensions(1) = 1 107 108* **** reset np_dimensions(3) if larger than nbrill **** 109 nbrill = control_num_kvectors() 110 if (np_dimensions(3).gt.nbrill) np_dimensions(3) = nbrill 111 112* **** reset np_dimensions(3) if it is not a multiple of np **** 113 do while ((mod(np,np_dimensions(3)).ne.0).and. 114 > (np_dimensions(3).gt.1)) 115 np_dimensions(3) = np_dimensions(3) - 1 116 end do 117* **** reset np_dimensions(2) if it is not a multiple of np2 **** 118 np = np/np_dimensions(3) 119 do while ((mod(np,np_dimensions(2)).ne.0).and. 120 > (np_dimensions(2).gt.1)) 121 np_dimensions(2) = np_dimensions(2) - 1 122 end do 123 124 !*** temporary restriction until ne parallelized in band *** 125 if ((code.eq.5) .or. 126 > (code.eq.10).or. 127 > (code.eq.13).or. 128 > (code.eq.14)) np_dimensions(2) = 1 129 130 np_dimensions(1) = np/(np_dimensions(2)*np_dimensions(3)) 131 132 endif 133#endif 134 135* **** get balance mapping **** 136 if (.not.btdb_get(rtdb,'nwpw:balance',mt_log,1,balance)) then 137 balance = .true. 138 end if 139 140* **** get parallel io **** 141 if (.not.btdb_get(rtdb,'nwpw:parallel_io',mt_log,1,pio)) then 142 pio = .false. 143 end if 144 !pio = pio.and.(np_dimensions(2).gt.1) 145 146* **** get fast_erf **** 147 if (.not.btdb_get(rtdb,'nwpw:fast_erf',mt_log,1,fast_erf)) then 148 fast_erf = .false. 149 end if 150 151* **** get fmm **** 152 if (.not.btdb_get(rtdb,'nwpw:fmm',mt_log,1,fmm)) then 153 fmm = .false. 154 end if 155 if (.not.btdb_get(rtdb,'nwpw:fmm_lmax',mt_int,1,fmm_lmax)) then 156 fmm_lmax = 10 157 end if 158 if (.not.btdb_get(rtdb,'nwpw:fmm_lr',mt_int,1,fmm_lr)) then 159 fmm_lr = 1 160 end if 161 162* **** get periodic_dipole **** 163 if (.not.btdb_get(rtdb,'nwpw:periodic_dipole', 164 > mt_log,1,periodic_dipole)) then 165 periodic_dipole = .false. 166 end if 167 168* **** get smooth_cuoff **** 169 smooth_cutoff = .true. 170 if (.not.btdb_get(rtdb,'nwpw:smooth_cutoff', 171 > mt_dbl,2,smooth_cutoff_values)) then 172 smooth_cutoff = .false. 173 smooth_cutoff_values(1) = 2.0d0 174 smooth_cutoff_values(2) = 4.0d0 175 end if 176 177* **** get hess_model mapping **** 178 if (.not.btdb_get(rtdb,'nwpw:hess_model',mt_log, 179 > 1,hess_model)) then 180 hess_model = .false. 181 end if 182 183 184 185* ********************************* 186* **** cpsd and band_sd: stuff **** 187* ********************************* 188 if ((code.eq.1).or.(code.eq.13)) then 189 if (.not.btdb_cget(rtdb,'cpsd:cell_name',1,cell_name)) then 190 cell_name = 'cell_default' 191 end if 192 193 if (.not.btdb_cget(rtdb,'cpsd:input_wavefunction_filename', 194 > 1,input_wavefunction_filename)) then 195 if (.not.btdb_cget(rtdb,'pspw:input vectors', 196 > 1,input_wavefunction_filename)) then 197 call util_file_prefix('movecs', 198 > input_wavefunction_filename) 199 end if 200 end if 201 202 if (.not.btdb_cget(rtdb,'cpsd:output_wavefunction_filename', 203 > 1,output_wavefunction_filename)) then 204 if (.not.btdb_cget(rtdb,'pspw:output vectors', 205 > 1,output_wavefunction_filename)) then 206 call util_file_prefix('movecs', 207 > output_wavefunction_filename) 208 end if 209 end if 210 211 212 if (.not.btdb_cget(rtdb,'cpsd:exchange_correlation', 213 > 1,exchange_correlation)) 214 > exchange_correlation = 'vosko' 215 216!$OMP single 217 if (nwpwxc_rtdb_load(rtdb,"dft")) then 218c call nwpwxc_print() 219 has_disp = nwpwxc_has_disp() 220 endif 221!$OMP end single 222 223#include "control_gga.fh" 224 225 if(.not.btdb_get(rtdb,'cpsd:geometry_optimize',mt_log,1,move)) 226 > move = .false. 227 if(.not.btdb_get(rtdb,'cpsd:fractional_coordinates', 228 > mt_log,1,frac_coord)) 229 > frac_coord = .false. 230 if(.not.btdb_get(rtdb,'cpsd:npsp',mt_int,1,npsp)) 231 > npsp=0 232 if(.not.btdb_get(rtdb,'cpsd:fake_mass',mt_dbl,1,fake_mass)) 233 > fake_mass = 400000.0d0 234 if(.not.btdb_get(rtdb,'cpsd:time_step',mt_dbl,1,time_step)) 235 > time_step = 5.8d0 236 if(.not.btdb_get(rtdb,'cpsd:loop',mt_int,2,loop)) then 237 loop(1)=10 238 loop(2)=100 239 end if 240 if(.not.btdb_get(rtdb,'cpsd:tolerances',mt_dbl,3,tolerances)) 241 >then 242 tolerances(1) = 1.0d-9 243 tolerances(2) = 1.0d-9 244 tolerances(3) = 1.0d-4 245 end if 246 scaling(1) = 0.0d0 247 scaling(2) = 0.0d0 248 249 if(.not.btdb_get(rtdb,'cpsd:ecut',mt_dbl,1,ecut)) 250 > ecut=9000.0d0 251 if(.not.btdb_get(rtdb,'cpsd:wcut',mt_dbl,1,wcut)) 252 > wcut=ecut 253 if(.not.btdb_get(rtdb,'cpsd:rcut',mt_dbl,1,rcut)) 254 > rcut = 0.0d0 255 if(.not.btdb_get(rtdb,'cpsd:ncut',mt_int,1,ncut)) 256 > ncut = 1 257 if(.not.btdb_get(rtdb,'cpsd:mult',mt_int,1,multiplicity)) 258 > multiplicity = 1 259 if(.not.btdb_get(rtdb,'cpsd:ispin',mt_int,1,ispin)) 260 > ispin=1 261 262 value3 = .false. 263 if (.not.btdb_get(rtdb,'nwpw:dof_rotation',mt_log,1,dof_rotation)) 264 > then 265 dof_rotation = .false. 266 value3= .true. 267 end if 268 269 if (.not.btdb_get(rtdb,'nwpw:rotation',mt_log,1,rotation)) then 270 rotation = .true. 271 else 272 if (value3) dof_rotation = rotation 273 end if 274 275 276 if (.not.btdb_get(rtdb,'nwpw:spin_orbit',mt_log,1,spin_orbit)) 277 > spin_orbit=.false. 278 if (spin_orbit) ispin=2 279 280 281* **************************************** 282* **** cpmd and band_cpmd code: stuff **** 283* **************************************** 284 else if ((code.eq.2) .or. (code.eq.14)) then 285 if (.not.btdb_cget(rtdb,'cpmd:cell_name',1,cell_name)) then 286 cell_name = 'cell_default' 287 end if 288 289 if(.not.btdb_cget(rtdb,'cpmd:input_wavefunction_filename', 290 > 1,input_wavefunction_filename)) then 291 if(.not.btdb_cget(rtdb,'pspw:input vectors', 292 > 1,input_wavefunction_filename)) then 293 call util_file_prefix('movecs', 294 > input_wavefunction_filename) 295 end if 296 end if 297 298 if(.not.btdb_cget(rtdb,'cpmd:output_wavefunction_filename', 299 > 1,output_wavefunction_filename)) then 300 if(.not.btdb_cget(rtdb,'pspw:output vectors', 301 > 1,output_wavefunction_filename)) then 302 call util_file_prefix('movecs', 303 > output_wavefunction_filename) 304 end if 305 end if 306 307 if(.not.btdb_cget(rtdb,'cpmd:input_v_wavefunction_filename', 308 > 1,input_v_wavefunction_filename)) then 309 if(.not.btdb_cget(rtdb,'pspw:input vvectors', 310 > 1,input_v_wavefunction_filename)) then 311 call util_file_prefix('vmovecs', 312 > input_v_wavefunction_filename) 313 end if 314 end if 315 316 if(.not.btdb_cget(rtdb,'cpmd:output_v_wavefunction_filename', 317 > 1,output_v_wavefunction_filename)) then 318 if(.not.btdb_cget(rtdb,'pspw:output vvectors', 319 > 1,output_v_wavefunction_filename)) then 320 call util_file_prefix('vmovecs', 321 > output_v_wavefunction_filename) 322 end if 323 end if 324 325 if(.not.btdb_cget(rtdb,'cpmd:xyz_filename', 326 > 1,xyz_filename)) 327 > call util_file_prefix('xyz',xyz_filename) 328! ind = index(exchange_correlation,' ') - 1 329 if(.not.btdb_cget(rtdb,'cpmd:exchange_correlation', 330 > 1,exchange_correlation)) 331 > exchange_correlation = 'vosko' 332 333!$OMP single 334 if (nwpwxc_rtdb_load(rtdb,"dft")) then 335c call nwpwxc_print() 336 has_disp = nwpwxc_has_disp() 337 endif 338!$OMP end single 339 340#include "control_gga.fh" 341 342 343 if(.not. btdb_get(rtdb,'cpmd:geometry_optimize',mt_log,1,move)) 344 > move = .true. 345 if(.not. btdb_get(rtdb,'cpmd:fractional_coordinates', 346 > mt_log,1,frac_coord)) 347 > frac_coord = .false. 348 if(.not. btdb_get(rtdb,'cpmd:npsp',mt_int,1,npsp)) 349 > npsp = 0 350 if(.not. btdb_get(rtdb,'cpmd:fake_mass',mt_dbl,1,fake_mass)) 351 > fake_mass = 800.0d0 352 if(.not. btdb_get(rtdb,'cpmd:time_step',mt_dbl,1,time_step)) 353 > time_step = 5.0d0 354 if(.not. btdb_get(rtdb,'cpmd:loop',mt_int,2,loop)) then 355 loop(1) = 10 356 loop(2) = 100 357 end if 358 if(.not. btdb_get(rtdb,'cpmd:scaling',mt_dbl,2,scaling)) then 359 scaling(1) = 1.0d0 360 scaling(2) = 1.0d0 361 end if 362 363 tolerances(1) = 0.0d0 364 tolerances(2) = 0.0d0 365 tolerances(3) = 0.0d0 366 if(.not. btdb_get(rtdb,'cpmd:ecut',mt_dbl,1,ecut)) 367 > ecut=9000.0d0 368 if(.not. btdb_get(rtdb,'cpmd:wcut',mt_dbl,1,wcut)) 369 > wcut = ecut 370 if(.not. btdb_get(rtdb,'cpmd:rcut',mt_dbl,1,rcut)) 371 > rcut = 0.0d0 372 if(.not. btdb_get(rtdb,'cpmd:ncut',mt_int,1,ncut)) 373 > ncut = 1 374 SA = .true. 375 if (.not.btdb_get(rtdb,'cpmd:sa_decay',mt_dbl,2,sa_decay)) then 376 SA = .false. 377 sa_decay(1) = 1.0d0 378 sa_decay(2) = 1.0d0 379 end if 380 381 if (.not.btdb_get(rtdb,'nwpw:dipole_motion',mt_log, 382 > 1,dipole_motion)) 383 > dipole_motion = .false. 384 385 if (.not.btdb_get(rtdb,'cpmd:fei',mt_log,1,fei)) 386 > fei = .false. 387 388 if (.not.btdb_get(rtdb,'cpmd:fei_quench',mt_log,1,fei_quench)) 389 > fei_quench = .false. 390 391 value3 = .false. 392 if (.not.btdb_get(rtdb,'nwpw:dof_rotation',mt_log,1,dof_rotation)) 393 > then 394 dof_rotation = .false. 395 value3= .true. 396 end if 397 398 if (.not.btdb_get(rtdb,'nwpw:rotation',mt_log,1,rotation)) then 399 rotation = .true. 400 else 401 if (value3) dof_rotation = rotation 402 end if 403 404 405* **** get thermostat information **** 406 if (.not.btdb_get(rtdb,'cpmd:nose',mt_log,1,nose)) 407 > nose = .false. 408 if (.not.btdb_get(rtdb,'cpmd:Pe',mt_dbl,1,Pe)) 409 > Pe = 1200.0d0 410 if (.not.btdb_get(rtdb,'cpmd:Te',mt_dbl,1,Te)) 411 > Te = 298.15d0 412 if (.not.btdb_get(rtdb,'cpmd:Pr',mt_dbl,1,Pr)) 413 > Pr = 1200.0d0 414 if (.not.btdb_get(rtdb,'cpmd:Tr',mt_dbl,1,Tr)) 415 > Tr = 298.15d0 416 417 if (.not.btdb_get(rtdb,'nwpw:spin_orbit',mt_log,1,spin_orbit)) 418 > spin_orbit=.false. 419 if (spin_orbit) ispin=2 420 421 422* ************************************ 423* **** cgsd: stuff or paw: stuff **** 424* ************************************ 425 else if ((code.eq.3).or.(code.eq.8).or. 426 > (code.eq.11).or.(code.eq.12).or.(code.eq.15)) then 427 if (.not.btdb_cget(rtdb,'cgsd:cell_name',1,cell_name)) then 428 cell_name = 'cell_default' 429 end if 430c 431c **** Figure input/output MO vectors **** 432c 433 if (.not. btdb_cget(rtdb, 'pspw:input vectors', 434 > 1,input_wavefunction_filename)) then 435 input_wavefunction_filename = 'atomic' 436 end if 437 438 if (.not. btdb_cget(rtdb, 'pspw:output vectors', 439 > 1,output_wavefunction_filename)) then 440 output_wavefunction_filename = ' ' 441 442 if (output_wavefunction_filename.eq.' ')then 443 if (input_wavefunction_filename.eq.'atomic')then 444 call util_file_prefix('movecs', 445 > output_wavefunction_filename) 446 else 447 output_wavefunction_filename=input_wavefunction_filename 448 endif 449 endif 450 if (input_wavefunction_filename.eq.'atomic')then 451 input_wavefunction_filename = output_wavefunction_filename 452 end if 453 end if 454 455 !**** determing input_ewavefunction_filename **** 456 if (.not.btdb_cget(rtdb, 'pspw:input evectors', 457 > 1,input_ewavefunction_filename)) then 458 if (.not.btdb_cget(rtdb,'cgsd:input_ewavefunction_filename', 459 > 1,input_ewavefunction_filename)) 460 > call util_file_prefix('emovecs',input_ewavefunction_filename) 461 end if 462 463 464 !**** determing output_ewavefunction_filename **** 465 if (.not.btdb_cget(rtdb,'pspw:output evectors', 466 > 1,output_ewavefunction_filename)) then 467 if (.not.btdb_cget(rtdb,'cgsd:output_ewavefunction_filename', 468 > 1,output_ewavefunction_filename)) then 469 output_ewavefunction_filename=input_ewavefunction_filename 470 end if 471 end if 472 473 474 if (code.eq.11) then 475 476 if(.not.btdb_cget(rtdb,'pspw:input vvectors', 477 > 1,input_v_wavefunction_filename)) then 478 if(.not.btdb_cget(rtdb,'cpmd:input_v_wavefunction_filename', 479 > 1,input_v_wavefunction_filename)) then 480 call util_file_prefix('vmovecs', 481 > input_v_wavefunction_filename) 482 end if 483 end if 484 485 if(.not.btdb_cget(rtdb,'pspw:output vvectors', 486 > 1,output_v_wavefunction_filename)) then 487 if(.not.btdb_cget(rtdb,'cpmd:output_v_wavefunction_filename', 488 > 1,output_v_wavefunction_filename)) then 489 490 call util_file_prefix('vmovecs', 491 > output_v_wavefunction_filename) 492 end if 493 end if 494 end if 495 496 497 if (.not.btdb_cget(rtdb,'cgsd:exchange_correlation', 498 > 1,exchange_correlation)) 499 > exchange_correlation = 'vosko' 500 501!$OMP single 502 if (nwpwxc_rtdb_load(rtdb,"dft")) then 503c call nwpwxc_print() 504 has_disp = nwpwxc_has_disp() 505 endif 506!$OMP end single 507 508#include "control_gga.fh" 509 510* ***** motion filenames **** 511 if(.not.btdb_cget(rtdb,'nwpw:xyz_filename', 512 > 1,xyz_filename)) 513 > call util_file_prefix('xyz',xyz_filename) 514 515 516* **** set Kohn-Sham scf parameters *** 517 if (.not. btdb_get(rtdb,'nwpw:ks_alpha',mt_dbl,1,ks_alpha)) 518 > ks_alpha = 0.25d0 519 if (.not. btdb_get(rtdb,'nwpw:fractional_alpha', 520 > mt_dbl,1,fractional_alpha)) 521 > fractional_alpha = 1.00d0 522 if (.not.btdb_get(rtdb,'nwpw:scf_algorithm', 523 > mt_int,1,scf_algorithm)) 524 > scf_algorithm = 3 525 526 if (.not.btdb_get(rtdb,'nwpw:diis_histories', 527 > mt_int,1,diis_histories)) 528 > diis_histories = 15 529 530 if (.not.btdb_get(rtdb,'nwpw:ks_algorithm', 531 > mt_int,1,ks_algorithm)) 532 > ks_algorithm = 0 533 if (.not.btdb_get(rtdb,'nwpw:kerker_g0', 534 > mt_dbl,1,kerker_g0)) 535 > kerker_g0 = 0.0d0 536 if (.not.btdb_get(rtdb,'nwpw:precondition', 537 > mt_log,1,precondition)) 538 > precondition = .false. 539 540* **** set maxit_orb maxit_orbs *** 541 if (.not.btdb_get(rtdb, 542 > 'nwpw:ks_maxit_orb',mt_int,1,maxit_orb)) 543 > maxit_orb = 5 544 if (.not.btdb_get(rtdb, 545 > 'nwpw:ks_maxit_orbs',mt_int,1,maxit_orbs)) 546 > maxit_orbs = 1 547 548 549 if (.not.btdb_get(rtdb,'cgsd:npsp',mt_int,1,npsp)) 550 > npsp = 0 551 if (.not.btdb_get(rtdb,'cgsd:fake_mass',mt_dbl,1,fake_mass)) 552 > fake_mass = 400000.0d0 553 if (.not.btdb_get(rtdb,'cgsd:time_step',mt_dbl,1,time_step)) 554 > time_step = 5.8d0 555 if (.not.btdb_get(rtdb,'cgsd:loop',mt_int,2,loop)) then 556 loop(1) = 10 557 loop(2) = 100 558 end if 559 if (.not.btdb_get(rtdb,'cgsd:tolerances',mt_dbl,3,tolerances))then 560 tolerances(1) = 1.0d-7 561 tolerances(2) = 1.0d-7 562 tolerances(3) = 1.0d-4 563 end if 564 565 if(.not. btdb_get(rtdb,'nwpw:scaling',mt_dbl,2,scaling)) then 566 scaling(1) = 1.0d0 567 scaling(2) = 1.0d0 568 end if 569 570 if (.not.btdb_get(rtdb,'cgsd:fractional_coordinates', 571 > mt_log,1,frac_coord)) 572 > frac_coord = .false. 573 574 if (.not.btdb_get(rtdb,'cgsd:ecut',mt_dbl,1,ecut)) 575 > ecut=9000.0d0 576 if (.not.btdb_get(rtdb,'cgsd:wcut',mt_dbl,1,wcut)) 577 > wcut = ecut 578 if (.not.btdb_get(rtdb,'cgsd:rcut',mt_dbl,1,rcut)) 579 > rcut = 0.0d0 580 if (.not.btdb_get(rtdb,'cgsd:ncut',mt_int,1,ncut)) 581 > ncut = 1 582 if (.not.btdb_get(rtdb,'cgsd:mult',mt_int,1,multiplicity)) 583 > multiplicity = 1 584 if (.not.btdb_get(rtdb,'cgsd:ispin',mt_int,1,ispin)) 585 > ispin = 1 586 587 588* **** BO parameterss *** 589 if (.not.btdb_get(rtdb,'nwpw:bo_steps',mt_int,2,bo_steps)) then 590 bo_steps(1) = 10 591 bo_steps(2) = 100 592 end if 593 if (.not.btdb_get(rtdb,'nwpw:bo_time_step',mt_dbl,1,bo_time_step)) 594 > bo_time_step = time_step 595 if (.not.btdb_get(rtdb,'nwpw:bo_algorithm',mt_int,1,bo_algorithm)) 596 > bo_algorithm = 0 597 if (.not.btdb_get(rtdb,'nwpw:bo_fake_mass',mt_dbl,1,bo_fake_mass)) 598 > bo_fake_mass = 500.0d0 599 600 value3 = .false. 601 if (.not.btdb_get(rtdb,'nwpw:dof_rotation',mt_log,1,dof_rotation)) 602 > then 603 dof_rotation = .false. 604 value3= .true. 605 end if 606 607 if (.not.btdb_get(rtdb,'nwpw:rotation',mt_log,1,rotation)) then 608 rotation = .true. 609 else 610 if (value3) dof_rotation = rotation 611 end if 612 613 614 SA = .true. 615 if (.not.btdb_get(rtdb,'nwpw:sa_decay',mt_dbl,2,sa_decay)) then 616 SA = .false. 617 sa_decay(1) = 1.0d0 618 sa_decay(2) = 1.0d0 619 end if 620 621 if (.not.btdb_get(rtdb,'nwpw:dipole_motion',mt_log, 622 > 1,dipole_motion)) 623 > dipole_motion = .false. 624 625 if (.not.btdb_get(rtdb,'nwpw:fei',mt_log,1,fei)) 626 > fei = .false. 627 628 if (.not.btdb_get(rtdb,'nwpw:fei_quench',mt_log,1,fei_quench)) 629 > fei_quench = .false. 630 631 632 633* **** get thermostat information **** 634 if (.not.btdb_get(rtdb,'nwpw:nose',mt_log,1,nose)) 635 > nose = .false. 636 if (.not.btdb_get(rtdb,'nwpw:Pe',mt_dbl,1,Pe)) 637 > Pe = 1200.0d0 638 if (.not.btdb_get(rtdb,'nwpw:Te',mt_dbl,1,Te)) 639 > Te = 298.15d0 640 if (.not.btdb_get(rtdb,'nwpw:Pr',mt_dbl,1,Pr)) 641 > Pr = 1200.0d0 642 if (.not.btdb_get(rtdb,'nwpw:Tr',mt_dbl,1,Tr)) 643 > Tr = 298.15d0 644 645 646* *************************** 647* **** pspw_dplot: stuff **** 648* *************************** 649 else if (code.eq.4) then 650 if (.not.btdb_cget(rtdb,'cgsd:cell_name',1,cell_name)) then 651 cell_name = 'cell_default' 652 end if 653 value = .true. 654 if (.not.btdb_cget(rtdb,'pspw_dplot:wavefunction_filename', 655 > 1,input_wavefunction_filename)) 656 > call util_file_prefix('movecs',input_wavefunction_filename) 657 call psi_get_header(i,ngrid,unita,ispin0,ne) 658 call dcopy(9,unita,1,unita_frozen,1) 659 if (i.eq.3) boundry = 'periodic' 660 if (i.eq.4) boundry = 'aperiodic' 661 662* **** dummy variables **** 663 move = .false. 664 frac_coord = .false. 665 gga = 0 666 fake_mass = 400000.0d0 667 time_step = 5.8d0 668 loop(1) = 0 669 loop(2) = 0 670 tolerances(1) = 1.0d-9 671 tolerances(2) = 1.0d-9 672 tolerances(3) = 1.0d-4 673 if(.not.btdb_get(rtdb,'cgsd:ecut',mt_dbl,1,ecut)) ecut=9000.0d0 674 if(.not.btdb_get(rtdb,'cgsd:wcut',mt_dbl,1,wcut)) wcut=ecut 675 rcut = 0.0d0 676 ncut = 1 677 npsp = 0 678 679c control_read = value 680c return 681 682* ********************* 683* **** band: stuff **** 684* ********************* 685 else if (code.eq.5) then 686 687 if (.not.btdb_cget(rtdb,'band:cell_name',1,cell_name)) then 688 cell_name = 'cell_default' 689 end if 690 691c **** Figure input/output MO vectors **** 692c 693 if (.not.btdb_cget(rtdb,'pspw:input vectors', 694 > 1,input_wavefunction_filename)) then 695 input_wavefunction_filename = 'atomic' 696 end if 697 698 if (.not.btdb_cget(rtdb, 'pspw:output vectors', 699 > 1,output_wavefunction_filename)) then 700 output_wavefunction_filename = ' ' 701 702 if (output_wavefunction_filename.eq.' ')then 703 if (input_wavefunction_filename.eq.'atomic') then 704 call util_file_prefix('movecs', 705 > output_wavefunction_filename) 706 else 707 output_wavefunction_filename=input_wavefunction_filename 708 endif 709 end if 710 end if 711 712 if (input_wavefunction_filename.eq.'atomic')then 713 input_wavefunction_filename = output_wavefunction_filename 714 end if 715 716 if (.not.btdb_cget(rtdb,'cgsd:input_ewavefunction_filename', 717 > 1,input_ewavefunction_filename)) 718 > call util_file_prefix('emovecs',input_ewavefunction_filename) 719 720 if (.not.btdb_cget(rtdb,'cgsd:output_ewavefunction_filename', 721 > 1,output_ewavefunction_filename)) 722 > call util_file_prefix('emovecs',output_ewavefunction_filename) 723 724 725 if (.not.btdb_cget(rtdb,'band:exchange_correlation', 726 > 1,exchange_correlation)) 727 > exchange_correlation = 'vosko' 728 729!$OMP single 730 if (nwpwxc_rtdb_load(rtdb,"dft")) then 731c call nwpwxc_print() 732 has_disp = nwpwxc_has_disp() 733 endif 734!$OMP end single 735 736#include "control_gga.fh" 737 738 739 if(.not.btdb_get(rtdb,'band:geometry_optimize',mt_log,1,move)) 740 > move = .false. 741 if (.not.btdb_get(rtdb,'band:fake_mass',mt_dbl,1,fake_mass)) 742 > fake_mass = 400000.0d0 743 if (.not.btdb_get(rtdb,'band:time_step',mt_dbl,1,time_step)) 744 > time_step = 5.8d0 745 if (.not.btdb_get(rtdb,'band:loop',mt_int,2,loop)) then 746 loop(1) = 10 747 loop(2) = 100 748 end if 749 if(.not.btdb_get(rtdb,'band:tolerances',mt_dbl,3,tolerances)) then 750 tolerances(1) = 1.0d-7 751 tolerances(2) = 1.0d-7 752 tolerances(3) = 1.0d-4 753 end if 754 scaling(1) = 0.0d0 755 scaling(2) = 0.0d0 756 757 if (.not.btdb_get(rtdb,'band:ecut',mt_dbl,1,ecut)) 758 > ecut = 9000.0d0 759 if (.not.btdb_get(rtdb,'band:wcut',mt_dbl,1,wcut)) 760 > wcut = ecut 761 if (.not.btdb_get(rtdb,'band:rcut',mt_dbl,1,rcut)) 762 > rcut = 0.0d0 763 if (.not.btdb_get(rtdb,'band:ncut',mt_int,1,ncut)) 764 > ncut = 1 765 if (.not.btdb_get(rtdb,'band:mult',mt_int,1,multiplicity)) 766 > multiplicity = 1 767 if (.not.btdb_get(rtdb,'band:ispin',mt_int,1,ispin)) 768 > ispin = 1 769 770 if (.not.btdb_get(rtdb,'band:spin_orbit',mt_log,1,spin_orbit)) 771 > spin_orbit=.false. 772 if (spin_orbit) ispin=2 773 774 775 776* **** set Kohn-Sham scf parameters *** 777 if (.not. btdb_get(rtdb,'nwpw:ks_alpha',mt_dbl,1,ks_alpha)) 778 > ks_alpha = 0.25d0 779 if (.not. btdb_get(rtdb,'nwpw:fractional_alpha', 780 > mt_dbl,1,fractional_alpha)) 781 > fractional_alpha = 1.00d0 782 if (.not.btdb_get(rtdb,'nwpw:scf_algorithm', 783 > mt_int,1,scf_algorithm)) 784 > scf_algorithm = 3 785 786 if (.not.btdb_get(rtdb,'nwpw:diis_histories', 787 > mt_int,1,diis_histories)) 788 > diis_histories = 15 789 790 if (.not.btdb_get(rtdb,'nwpw:ks_algorithm', 791 > mt_int,1,ks_algorithm)) 792 > ks_algorithm = 0 793 if (.not.btdb_get(rtdb,'nwpw:kerker_g0', 794 > mt_dbl,1,kerker_g0)) 795 > kerker_g0 = 0.0d0 796 if (.not.btdb_get(rtdb,'nwpw:precondition', 797 > mt_log,1,precondition)) 798 > precondition = .false. 799 800 801* **** set maxit_orb maxit_orbs *** 802 if (.not.btdb_get(rtdb, 803 > 'nwpw:ks_maxit_orb',mt_int,1,maxit_orb)) 804 > maxit_orb = 5 805 if (.not.btdb_get(rtdb, 806 > 'nwpw:ks_maxit_orbs',mt_int,1,maxit_orbs)) 807 > maxit_orbs = 1 808 809 810* *********************** 811* **** paw_sd: stuff **** 812* *********************** 813 else if (code.eq.6) then 814 815 if (.not.btdb_get(rtdb,'cgsd:geometry_optimize',mt_log,1,move)) 816 > move = .false. 817 818 if (.not.btdb_cget(rtdb,'cgsd:cell_name',1,cell_name)) then 819 cell_name = 'cell_default' 820 end if 821 if (.not.btdb_cget(rtdb,'cgsd:input_wavefunction_filename', 822 > 1,input_wavefunction_filename)) 823 > call util_file_prefix('movecs',input_wavefunction_filename) 824 if (.not.btdb_cget(rtdb,'cgsd:output_wavefunction_filename', 825 > 1,output_wavefunction_filename)) 826 > call util_file_prefix('movecs',output_wavefunction_filename) 827 if (.not.btdb_cget(rtdb,'cgsd:exchange_correlation', 828 > 1,exchange_correlation)) 829 > exchange_correlation = 'vosko' 830 831 if (inp_compare(.false.,exchange_correlation,'vosko')) then 832 gga = 0 833 else if (inp_compare(.false.,exchange_correlation,'lda')) then 834 gga = 0 835 else if (inp_compare(.false.,exchange_correlation,'svwn5')) then 836 gga = 0 837 else if (inp_compare(.false.,exchange_correlation,'pbe96')) then 838 gga = 10 839 else if (inp_compare(.false.,exchange_correlation,'blyp')) then 840 gga = 11 841 else if (inp_compare(.false.,exchange_correlation,'revpbe')) then 842 gga = 12 843 else 844 gga = 0 845 end if 846 847!$OMP single 848 if (nwpwxc_rtdb_load(rtdb,"dft")) then 849c call nwpwxc_print() 850 has_disp = nwpwxc_has_disp() 851 endif 852!$OMP end single 853 854 if (.not.btdb_get(rtdb,'cgsd:npsp',mt_int,1,npsp)) 855 > npsp = 0 856 if (.not.btdb_get(rtdb,'cgsd:fake_mass',mt_dbl,1,fake_mass)) 857 > fake_mass = 400000.0d0 858 if (.not.btdb_get(rtdb,'cgsd:time_step',mt_dbl,1,time_step)) 859 > time_step = 5.8d0 860 if (.not.btdb_get(rtdb,'cgsd:loop',mt_int,2,loop)) then 861 loop(1) = 10 862 loop(2) = 100 863 end if 864 if(.not.btdb_get(rtdb,'cgsd:tolerances',mt_dbl,3,tolerances)) then 865 tolerances(1) = 1.0d-7 866 tolerances(2) = 1.0d-7 867 tolerances(3) = 1.0d-4 868 end if 869 scaling(1) = 0.0d0 870 scaling(2) = 0.0d0 871 if (.not.btdb_get(rtdb,'cgsd:ecut',mt_dbl,1,ecut)) 872 > ecut = 9000.0d0 873 if (.not.btdb_get(rtdb,'cgsd:wcut',mt_dbl,1,wcut)) 874 > wcut = ecut 875 if (.not.btdb_get(rtdb,'cgsd:rcut',mt_dbl,1,rcut)) 876 > rcut = 0.0d0 877 if (.not.btdb_get(rtdb,'cgsd:ncut',mt_int,1,ncut)) 878 > ncut = 1 879 if (.not.btdb_get(rtdb,'cgsd:mult',mt_int,1,multiplicity)) 880 > multiplicity = 1 881 if (.not.btdb_get(rtdb,'cgsd:ispin',mt_int,1,ispin)) 882 > ispin = 1 883 884 885* ************************* 886* **** paw_cpmd: stuff **** 887* ************************* 888 else if (code.eq.7) then 889 890 if (.not.btdb_cget(rtdb,'cpmd:cell_name',1,cell_name)) then 891 cell_name = 'cell_default' 892 end if 893 894 if (.not.btdb_cget(rtdb,'cpmd:input_wavefunction_filename', 895 > 1,input_wavefunction_filename)) 896 > call util_file_prefix('movecs',input_wavefunction_filename) 897 if (.not.btdb_cget(rtdb,'cpmd:output_wavefunction_filename', 898 > 1,output_wavefunction_filename)) 899 > call util_file_prefix('movecs',output_wavefunction_filename) 900 if (.not.btdb_cget(rtdb,'cpmd:input_v_wavefunction_filename', 901 > 1,input_v_wavefunction_filename)) 902 > call util_file_prefix('vmovecs',input_v_wavefunction_filename) 903 if (.not.btdb_cget(rtdb,'cpmd:output_v_wavefunction_filename', 904 > 1,output_v_wavefunction_filename)) 905 > call util_file_prefix('vmovecs',output_v_wavefunction_filename) 906 907 908* ***** motion filenames **** 909 if(.not.btdb_cget(rtdb,'cpmd:xyz_filename', 910 > 1,xyz_filename)) 911 > call util_file_prefix('xyz',xyz_filename) 912 913 914 if (.not.btdb_cget(rtdb,'cpmd:exchange_correlation', 915 > 1,exchange_correlation)) 916 > exchange_correlation = 'vosko' 917 918 if (inp_compare(.false.,exchange_correlation,'vosko')) then 919 gga = 0 920 else if (inp_compare(.false.,exchange_correlation,'lda')) then 921 gga = 0 922 else if (inp_compare(.false.,exchange_correlation,'svwn5')) then 923 gga = 0 924 else if (inp_compare(.false.,exchange_correlation,'pbe96')) then 925 gga = 10 926 else if (inp_compare(.false.,exchange_correlation,'blyp')) then 927 gga = 11 928 else if (inp_compare(.false.,exchange_correlation,'revpbe')) then 929 gga = 12 930 else 931 gga = 0 932 end if 933 934!$OMP single 935 if (nwpwxc_rtdb_load(rtdb,"dft")) then 936c call nwpwxc_print() 937 has_disp = nwpwxc_has_disp() 938 endif 939!$OMP end single 940 941 942 if (.not.btdb_get(rtdb,'cpmd:geometry_optimize',mt_log,1,move)) 943 > move = .false. 944 if (.not.btdb_get(rtdb,'cpmd:fractional_coordinates', 945 > mt_log,1,frac_coord)) 946 > frac_coord = .false. 947 if (.not.btdb_get(rtdb,'cpmd:npsp',mt_int,1,npsp)) 948 > npsp = 0 949 if (.not.btdb_get(rtdb,'cpmd:fake_mass',mt_dbl,1,fake_mass)) 950 > fake_mass = 800.0d0 951 if (.not.btdb_get(rtdb,'cpmd:time_step',mt_dbl,1,time_step)) 952 > time_step = 5.0d0 953 if (.not.btdb_get(rtdb,'cpmd:loop',mt_int,2,loop)) then 954 loop(1) = 10 955 loop(2) = 100 956 end if 957 if (.not.btdb_get(rtdb,'cpmd:scaling',mt_dbl,2,scaling)) then 958 scaling(1) = 1.0d0 959 scaling(2) = 1.0d0 960 end if 961 tolerances(1) = 0.0d0 962 tolerances(2) = 0.0d0 963 tolerances(3) = 0.0d0 964 if (.not.btdb_get(rtdb,'cpmd:ecut',mt_dbl,1,ecut)) 965 > ecut = 9000.0d0 966 if (.not.btdb_get(rtdb,'cpmd:wcut',mt_dbl,1,wcut)) 967 > wcut = ecut 968 if (.not.btdb_get(rtdb,'cpmd:rcut',mt_dbl,1,rcut)) 969 > rcut = 0.0d0 970 if (.not.btdb_get(rtdb,'cpmd:ncut',mt_int,1,ncut)) 971 > ncut = 1 972 973 SA = .true. 974 if (.not.btdb_get(rtdb,'cpmd:sa_decay',mt_dbl,2,sa_decay)) then 975 SA = .false. 976 sa_decay(1) = 1.0d0 977 sa_decay(2) = 1.0d0 978 end if 979 980 if (.not.btdb_get(rtdb,'nwpw:dipole_motion',mt_log, 981 > 1,dipole_motion)) 982 > dipole_motion = .false. 983 984 if (.not.btdb_get(rtdb,'cpmd:fei',mt_log,1,fei)) 985 > fei = .false. 986 987 if (.not.btdb_get(rtdb,'cpmd:fei_quench',mt_log,1,fei_quench)) 988 > fei_quench = .false. 989 990 value3 = .false. 991 if (.not.btdb_get(rtdb,'nwpw:dof_rotation',mt_log,1,dof_rotation)) 992 > then 993 dof_rotation = .false. 994 value3= .true. 995 end if 996 997 if (.not.btdb_get(rtdb,'nwpw:rotation',mt_log,1,rotation)) then 998 rotation = .true. 999 else 1000 if (value3) dof_rotation = rotation 1001 end if 1002 1003* **** get thermostat information **** 1004 if (.not.btdb_get(rtdb,'cpmd:nose',mt_log,1,nose)) 1005 > nose = .false. 1006 if (.not.btdb_get(rtdb,'cpmd:Pe',mt_dbl,1,Pe)) 1007 > Pe = 1200.0d0 1008 if (.not.btdb_get(rtdb,'cpmd:Te',mt_dbl,1,Te)) 1009 > Te = 298.15d0 1010 if (.not.btdb_get(rtdb,'cpmd:Pr',mt_dbl,1,Pr)) 1011 > Pr = 1200.0d0 1012 if (.not.btdb_get(rtdb,'cpmd:Tr',mt_dbl,1,Tr)) 1013 > Tr = 298.15d0 1014 1015 1016* ********************** 1017* **** pspw_wannier **** 1018* ********************** 1019 else if (code.eq.9) then 1020 1021 if (.not.btdb_cget(rtdb,'cgsd:cell_name',1,cell_name)) then 1022 cell_name = 'cell_default' 1023 end if 1024 1025c **** Figure input/output MO vectors **** 1026 1027 if (.not.btdb_cget(rtdb,'wannier:input vectors', 1028 > 1,input_wavefunction_filename)) then 1029 if (.not.btdb_cget(rtdb, 'pspw:input vectors', 1030 > 1,input_wavefunction_filename)) then 1031 input_wavefunction_filename = 'atomic' 1032 end if 1033 end if 1034 1035 if (.not.btdb_cget(rtdb,'wannier:output vectors', 1036 > 1,output_wavefunction_filename)) then 1037 if (.not.btdb_cget(rtdb,'pspw:output vectors', 1038 > 1,output_wavefunction_filename)) then 1039 output_wavefunction_filename = ' ' 1040 end if 1041 end if 1042 1043 if (output_wavefunction_filename.eq.' ')then 1044 if (input_wavefunction_filename.eq.'atomic')then 1045 call util_file_prefix('movecs',output_wavefunction_filename) 1046 else 1047 output_wavefunction_filename = input_wavefunction_filename 1048 endif 1049 endif 1050 1051 if (input_wavefunction_filename.eq.'atomic')then 1052 input_wavefunction_filename = output_wavefunction_filename 1053 end if 1054 1055 call psi_get_header(i,ngrid,unita,ispin0,ne) 1056 call dcopy(9,unita,1,unita_frozen,1) 1057 if (i.eq.3) boundry = 'periodic' 1058 if (i.eq.4) boundry = 'aperiodic' 1059 1060* **** dummy variables **** 1061 move = .false. 1062 frac_coord = .false. 1063 gga = 0 1064 fake_mass = 400000.0d0 1065 time_step = 5.8d0 1066 loop(1) = 0 1067 loop(2) = 0 1068 tolerances(1) = 1.0d-9 1069 tolerances(2) = 1.0d-9 1070 tolerances(3) = 1.0d-4 1071 if(.not.btdb_get(rtdb,'cgsd:ecut',mt_dbl,1,ecut)) ecut=9000.0d0 1072 if(.not.btdb_get(rtdb,'cgsd:wcut',mt_dbl,1,wcut)) wcut=ecut 1073 rcut = 0.0d0 1074 ncut = 0 1075 npsp = 0 1076 1077c control_read = value 1078c return 1079 1080* *************************** 1081* **** band_dplot: stuff **** 1082* *************************** 1083 else if (code.eq.10) then 1084 1085 if (.not.btdb_cget(rtdb,'band:cell_name',1,cell_name)) then 1086 cell_name = 'cell_default' 1087 end if 1088 1089 value = .true. 1090 if (.not.btdb_cget(rtdb,'band_dplot:wavefunction_filename', 1091 > 1,input_wavefunction_filename)) 1092 > call util_file_prefix('movecs',input_wavefunction_filename) 1093 call cpsi_get_header(i,ngrid,unita,ispin0,ne,nbrill) 1094 call dcopy(9,unita,1,unita_frozen,1) 1095 if (i.eq.5) boundry = 'periodic' 1096 1097* **** dummy variables **** 1098 move = .false. 1099 frac_coord = .false. 1100 gga = 0 1101 fake_mass = 400000.0d0 1102 time_step = 5.8d0 1103 loop(1) = 0 1104 loop(2) = 0 1105 tolerances(1) = 1.0d-9 1106 tolerances(2) = 1.0d-9 1107 tolerances(3) = 1.0d-4 1108 if(.not.btdb_get(rtdb,'band:ecut',mt_dbl,1,ecut)) ecut=9000.0d0 1109 if(.not.btdb_get(rtdb,'band:wcut',mt_dbl,1,wcut)) wcut=ecut 1110 rcut = 0.0d0 1111 ncut = 0 1112 npsp = 0 1113 1114c control_read = value 1115c return 1116 1117* ******************************** 1118* **** unknown but dont crash **** 1119* ******************************** 1120 else if (code.gt.99) then 1121 call nwpw_timing_end(50) 1122 control_read = value 1123 return 1124 1125* *************************** 1126* **** unknown code type **** 1127* *************************** 1128 else 1129 value = .false. 1130 write(*,*) "control_read: unknown code type:",code 1131 call nwpw_timing_end(50) 1132 control_read = value 1133 return 1134 end if 1135 1136* ***************************** 1137* ***** symmetry variables **** 1138* ***************************** 1139 if (.not.btdb_get(rtdb,'nwpw:symmetry',mt_int,1,symm_number)) 1140 > symm_number = 0 1141 1142* ********************** 1143* ***** cell: stuff **** 1144* ********************** 1145 l = index(cell_name,' ') - 1 1146 rtdb_unita = cell_name(1:l)//':unita' 1147 rtdb_unitaf = cell_name(1:l)//':unita_frozen' 1148 rtdb_ngrid = cell_name(1:l)//':ngrid' 1149 rtdb_boundry = cell_name(1:l)//':boundry' 1150 rtdb_ngrid_small = cell_name(1:l)//':ngrid_small' 1151 1152 1153* **** define unita and boundary **** 1154 if (.not.btdb_get(rtdb,rtdb_unita,mt_dbl,9,unita)) then 1155 call dcopy(9,0.0d0,0,unita,1) 1156 end if 1157 1158 if (.not.btdb_cget(rtdb,rtdb_boundry,1,boundry)) then 1159 boundry = 'periodic' 1160 end if 1161 call check_unita_for_default(rtdb,unita,rtdb_unita,cell_name) 1162 1163 1164* **** define unita_frozen **** 1165 if (.not.btdb_get(rtdb,'nwpw:frozen_lattice',mt_log,1,frozen)) 1166 > frozen = .true. 1167 1168 if (frozen) then 1169 if (.not.btdb_get(rtdb,rtdb_unitaf,mt_dbl,9,unita_frozen)) then 1170 call dcopy(9,unita,1,unita_frozen,1) 1171 value2 = btdb_parallel(.false.) 1172 if (taskid.eq.MASTER) then 1173 value = value.and. 1174 > btdb_put(rtdb,rtdb_unitaf,mt_dbl,9,unita_frozen) 1175 end if 1176 value2 = btdb_parallel(.true.) 1177 else 1178 if (.not.btdb_get(rtdb,'nwpw:frozen_lattice:thresh', 1179 > mt_dbl,1,thresh)) thresh = 0.05d0 1180 error = 0.0d0 1181 do l=1,3 1182 do i=1,3 1183 error = error + ((unita_frozen(i,l)-unita(i,l)))**2 1184 end do 1185 end do 1186 error = dsqrt(error) 1187 if (error.gt.thresh) then 1188 call dcopy(9,unita,1,unita_frozen,1) 1189 value2 = btdb_parallel(.false.) 1190 if (taskid.eq.MASTER) then 1191 value = value.and. 1192 > btdb_put(rtdb,rtdb_unitaf,mt_dbl,9,unita_frozen) 1193 end if 1194 value2 = btdb_parallel(.true.) 1195 end if 1196 end if 1197 else 1198 call dcopy(9,unita,1,unita_frozen,1) 1199 end if 1200 1201 1202* **** define ngrid based on unita and ecut **** 1203 if (.not.btdb_get(rtdb,rtdb_ngrid,mt_int,3,ngrid)) then 1204 if ((ecut.gt.5000.0d0).and.(wcut.gt.5000.0d0)) then 1205 call control_ecut_wcut_default(rtdb,ecut,wcut) 1206 else if ((ecut.gt.5000.0d0).and.(wcut.lt.5000.0d0)) then 1207 ecut = 2.0d0*wcut 1208 end if 1209 1210 !call control_ngrid_default(rtdb,unita,ecut,mapping,ngrid) 1211 call control_ngrid_default(rtdb,unita_frozen,ecut, 1212 > mapping,ngrid) 1213 !ngrid(1) = 32 1214 !ngrid(2) = 32 1215 !ngrid(3) = 32 1216 end if 1217 1218 1219* **** define ngrid_small **** 1220 has_ngrid_small = .false. 1221 ngrid_small(1) = 0 1222 ngrid_small(2) = 0 1223 ngrid_small(3) = 0 1224 if (btdb_get(rtdb,rtdb_ngrid_small,mt_int,3,ngrid_small)) 1225 > has_ngrid_small = .true. 1226 1227* *** set to false if wannier or dplot code *** 1228 if ((code.eq.4) .or.(code.eq.9).or.(code.eq.10)) 1229 > has_ngrid_small = .false. 1230 1231* **** set single_p;recision hfx **** 1232 if (.not.btdb_get(rtdb,'pspw:HFX_single_precision',mt_log,1, 1233 > single_precision_on)) 1234 > single_precision_on = .false. 1235 1236 1237* **** set fractional (smearing) parameters **** 1238 if (.not. 1239 > btdb_get(rtdb,'nwpw:fractional_orbitals',mt_int,2,frac_ne)) then 1240 frac_ne(1) = 0 1241 frac_ne(2) = 0 1242 end if 1243 fractional = (frac_ne(1).gt.0).or.(frac_ne(2).gt.0) 1244 if (.not.btdb_get(rtdb,'nwpw:fractional_temperature', 1245 > mt_dbl,1,frac_temperature)) then 1246 frac_temperature = 0.0d0 1247 end if 1248 if (.not.btdb_get(rtdb,'nwpw:fractional_smeartype', 1249 > mt_int,1,frac_smeartype)) then 1250 frac_smeartype = 0 1251 end if 1252 if (.not. btdb_get(rtdb,'nwpw:fractional_alpha', 1253 > mt_dbl,1,fractional_alpha)) then 1254 fractional_alpha = 1.00d0 1255 end if 1256 1257* **** set attenuation parameter **** 1258 if (.not.btdb_get(rtdb,'nwpw:attenuation',mt_dbl,1,attenuation)) 1259 > attenuation = 0.5d0 1260 1261* **** set preconditioning parameters Ep,Sp **** 1262 if (.not.btdb_get(rtdb,'nwpw:Eprecondition',mt_dbl,1,Ep)) 1263 > Ep = 20.0d0 1264 if (.not.btdb_get(rtdb,'nwpw:Sprecondition',mt_dbl,1,Sp)) 1265 > Sp = 200.0d0 1266 1267* **** set out of time variables **** 1268 est_step_time = -1 1269 est_finish_time = -1 1270 call current_second(cpu1_time) 1271 1272* **** set gram_schmidt *** 1273 gram_schmidt = .false. 1274 if (.not.btdb_get(rtdb, 1275 > 'nwpw:gram-schmidt',mt_log,1,gram_schmidt)) 1276 > gram_schmidt = .false. 1277 1278* **** set translation *** 1279 translation = control_allow_translation() 1280 dof_translation = translation 1281 1282* **** set two component pseudopotential 1283 two_comp_ppot = .false. 1284 1285 1286* **** ewald ngrid ************ 1287cccc set to some reasonable default ! 1288 if (.not.btdb_get(rtdb,'nwpw:ewald_ngrid',mt_int,3,ewald_grid)) 1289 > then 1290 ewald_grid(1)=ngrid(1) 1291 ewald_grid(2)=ngrid(2) 1292 ewald_grid(3)=ngrid(3) 1293 end if 1294 1295* **** reset mapping if slab and it doesnot fit **** 1296 if (mapping.eq.1) then 1297 if ((ngrid(2).ne.ngrid(3)).or. 1298 > (ngrid(3).lt.np_dimensions(1))) then 1299 mapping = 2 1300 end if 1301 end if 1302 1303 if (np_default) then 1304 if ( (code.eq.1).or. 1305 > (code.eq.2).or. 1306 > (code.eq.3).or. 1307 > (code.eq.9).or. 1308 > (code.eq.11)) then 1309 nmult(1) = 2 1310 nmult(2) = 3 1311 nmult(3) = 5 1312 do i=1,3 1313 do while ((np_dimensions(1).gt.ngrid(1)).and. 1314 > (mod(np_dimensions(1),nmult(i)).eq.0)) 1315 np_dimensions(1) = np_dimensions(1)/nmult(i) 1316 np_dimensions(2) = np_dimensions(2)*nmult(i) 1317 end do 1318 do while ((np_dimensions(1).gt.ngrid(2)).and. 1319 > (mod(np_dimensions(1),nmult(i)).eq.0)) 1320 np_dimensions(1) = np_dimensions(1)/nmult(i) 1321 np_dimensions(2) = np_dimensions(2)*nmult(i) 1322 end do 1323 do while ((np_dimensions(1).gt.ngrid(3)).and. 1324 > (mod(np_dimensions(1),nmult(i)).eq.0)) 1325 np_dimensions(1) = np_dimensions(1)/nmult(i) 1326 np_dimensions(2) = np_dimensions(2)*nmult(i) 1327 end do 1328 end do 1329 end if 1330 end if 1331 pio = pio.and.((np_dimensions(2).gt.1).or.(np_dimensions(3).gt.1)) 1332 1333* **** minimizer **** 1334 if (.not.btdb_get(rtdb,'nwpw:minimizer',mt_int,1,minimizer)) then 1335 minimizer = 1 ! make the default Grassmann cg 1336 end if 1337 1338 call nwpw_timing_end(50) 1339 control_read = value 1340 return 1341 end 1342* *********************************** 1343* * control_ewald_ngrid() 1344* *********************************** 1345 integer function control_ewald_ngrid(i) 1346 implicit none 1347#include "control.fh" 1348 integer i 1349 control_ewald_ngrid=ewald_grid(i) 1350 return 1351 end 1352* *********************************** 1353* * control_ewald_set_ngrid 1354* *********************************** 1355 subroutine ewald_set_ngrid(enx,eny,enz) 1356 implicit none 1357 integer enx,eny,enz 1358 1359#include "control.fh" 1360 1361 ewald_grid(1)=enx 1362 ewald_grid(2)=eny 1363 ewald_grid(3)=enz 1364 return 1365 end 1366 1367* *********************************** 1368* * * 1369* * control_ngrid_default * 1370* * * 1371* *********************************** 1372 subroutine control_ngrid_default(rtdb,unita,ecut,mapping,ngrid) 1373 implicit none 1374 integer rtdb 1375 real*8 unita(3,3),ecut 1376 integer mapping 1377 integer ngrid(3) 1378 1379#include "bafdecls.fh" 1380#include "btdb.fh" 1381#include "errquit.fh" 1382 1383* **** local variables **** 1384 real*8 unitg(3,3),omega 1385 real*8 gx,gy,gz 1386 real*8 xh,yh,zh 1387 integer control_set_ngrid 1388 external control_set_ngrid 1389 1390 call get_cube(unita,unitg,omega) 1391 1392 gx = unitg(1,1) 1393 gy = unitg(2,1) 1394 gz = unitg(3,1) 1395 xh = dsqrt(2.0d0*ecut/(gx*gx + gy*gy + gz*gz))+0.5d0 1396 1397 gx = unitg(1,2) 1398 gy = unitg(2,2) 1399 gz = unitg(3,2) 1400 yh = dsqrt(2.0d0*ecut/(gx*gx + gy*gy + gz*gz))+0.5d0 1401 1402 gx = unitg(1,3) 1403 gy = unitg(2,3) 1404 gz = unitg(3,3) 1405 zh = dsqrt(2.0d0*ecut/(gx*gx + gy*gy + gz*gz))+0.5d0 1406 1407 if (mapping.ge.2) then 1408c ngrid(1) = control_set_ngrid(2.0d0*xh,.true.) 1409c ngrid(2) = control_set_ngrid(2.0d0*yh,.false.) 1410c ngrid(3) = control_set_ngrid(2.0d0*zh,.false.) 1411 ngrid(1) = control_set_ngrid(2.0d0*xh,.true.) 1412 ngrid(2) = control_set_ngrid(2.0d0*yh,.true.) 1413 ngrid(3) = control_set_ngrid(2.0d0*zh,.true.) 1414 else 1415 ngrid(1) = control_set_ngrid(2.0d0*xh,.true.) 1416 ngrid(2) = control_set_ngrid(2.0d0*yh,.true.) 1417 ngrid(3) = control_set_ngrid(2.0d0*zh,.true.) 1418 if (ngrid(2).gt.ngrid(3)) then 1419 ngrid(3) = ngrid(2) 1420 else 1421 ngrid(2) = ngrid(3) 1422 end if 1423 end if 1424 1425 1426c* *** write unita to rtdb - should happen only once during a simulation *** 1427c if (.not.btdb_put(rtdb,rtdb_ngrid,mt_int,3,ngrid)) then 1428c call errquit('cannot write ngrid to rtdb',0,0) 1429c end if 1430 1431 return 1432 end 1433 1434* *********************************** 1435* * * 1436* * control_set_ngrid * 1437* * * 1438* *********************************** 1439* 1440* return n so that it is a multiple of 2,3,5,7 1441* 1442 integer function control_set_ngrid(x,mult2) 1443 implicit none 1444 real*8 x 1445 logical mult2 1446 1447* **** local variables **** 1448 integer nx,ntest 1449 integer nf2,nf3,nf5,nf7 1450 1451 integer factor_count2 1452 external factor_count2 1453 1454* **** find prime factors of nx *** 1455 nx = (x+0.5d0) !*** crude rounding 1456 if ((mult2).and.(mod(nx,2).ne.0)) nx = nx+1 1457 1458 nf2 = factor_count2(nx,2) 1459 nf3 = factor_count2(nx,3) 1460 nf5 = factor_count2(nx,5) 1461 nf7 = factor_count2(nx,7) 1462 ntest = (2**nf2) * (3**nf3) * (5**nf5) * (7**nf7) 1463 do while (nx .ne. ntest) 1464 nx = nx + 1 1465 if (mult2) nx = nx + 1 1466 nf2 = factor_count2(nx,2) 1467 nf3 = factor_count2(nx,3) 1468 nf5 = factor_count2(nx,5) 1469 nf7 = factor_count2(nx,7) 1470 ntest = (2**nf2) * (3**nf3) * (5**nf5) * (7**nf7) 1471 end do 1472 1473 control_set_ngrid = nx 1474 return 1475 end 1476 1477 integer function factor_count2(n,m) 1478 implicit none 1479 integer n,m 1480 integer f,nn 1481 1482 f = 0 1483 nn = n 1484 1485 do while (mod(nn,m).eq.0) 1486 nn = nn/m 1487 f = f + 1 1488 end do 1489 1490 factor_count2 = f 1491 return 1492 end 1493 1494 1495 1496 1497 1498* *********************************** 1499* * * 1500* * check_unita_for_default * 1501* * * 1502* *********************************** 1503 subroutine check_unita_for_default(rtdb,unita,rtdb_unita, 1504 > cell_name) 1505 implicit none 1506 integer rtdb 1507 real*8 unita(3,3) 1508 character*50 rtdb_unita 1509 character*50 cell_name 1510 1511 1512 1513#include "bafdecls.fh" 1514#include "btdb.fh" 1515#include "beom.fh" 1516#include "errquit.fh" 1517 1518* **** local variables **** 1519 integer MASTER,taskid 1520 parameter (MASTER=0) 1521 logical value,box_orient 1522 integer geom,box_type,l,isystype 1523 real*8 box_delta 1524 character*50 rtdb_name 1525 1526 value = (unita(1,1) .eq. 0.0d0).and. 1527 > (unita(2,1) .eq. 0.0d0).and. 1528 > (unita(3,1) .eq. 0.0d0).and. 1529 > (unita(1,2) .eq. 0.0d0).and. 1530 > (unita(2,2) .eq. 0.0d0).and. 1531 > (unita(3,2) .eq. 0.0d0).and. 1532 > (unita(1,3) .eq. 0.0d0).and. 1533 > (unita(2,3) .eq. 0.0d0).and. 1534 > (unita(3,3) .eq. 0.0d0) 1535 1536 if (value) then 1537 value = beom_create(geom,'geometry') 1538 value = value.and.beom_rtdb_load(rtdb,geom,'geometry') 1539 value = value.and.geom_amatrix_get(geom,unita) 1540 value = value.and.geom_systype_get(geom,isystype) 1541 if (.not. value) call errquit('cannot load geometry',0, 1542 & GEOM_ERR) 1543 1544 value = (unita(1,1) .eq. 1.0d0).and. 1545 > (unita(2,1) .eq. 0.0d0).and. 1546 > (unita(3,1) .eq. 0.0d0).and. 1547 > (unita(1,2) .eq. 0.0d0).and. 1548 > (unita(2,2) .eq. 1.0d0).and. 1549 > (unita(3,2) .eq. 0.0d0).and. 1550 > (unita(1,3) .eq. 0.0d0).and. 1551 > (unita(2,3) .eq. 0.0d0).and. 1552 > (unita(3,3) .eq. 1.0d0) 1553 if (value) then 1554 l = index(cell_name,' ') - 1 1555 rtdb_name = cell_name(1:l)//':box_delta' 1556 if (.not.btdb_get(rtdb,rtdb_name,mt_dbl,1,box_delta)) 1557 > box_delta = 5.0d0 1558 rtdb_name = cell_name(1:l)//':box_type' 1559 if (.not.btdb_get(rtdb,rtdb_name,mt_int,1,box_type)) 1560 > box_type = 1 1561 rtdb_name = cell_name(1:l)//':box_orient' 1562 if (.not.btdb_get(rtdb,rtdb_name,mt_log,1, 1563 > box_orient)) 1564 > box_orient = .true. 1565 call control_find_box(geom,box_type,box_orient,box_delta, 1566 > unita) 1567 1568* *** write unita to rtdb - should happen only once during a simulation *** 1569 call Parallel_taskid(taskid) 1570 value = btdb_parallel(.false.) 1571 if (taskid.eq.MASTER) then 1572 if (.not.btdb_put(rtdb,rtdb_unita,mt_dbl,9,unita)) then 1573 call errquit('cannot write unita to rtdb',0,0) 1574 end if 1575 end if 1576 value = btdb_parallel(.true.) 1577 1578 !unita(1,1) = 20.0d0 1579 !unita(2,1) = 0.0d0 1580 !unita(3,1) = 0.0d0 1581 !unita(1,2) = 0.0d0 1582 !unita(2,2) = 20.0d0 1583 !unita(3,2) = 0.0d0 1584 !unita(1,3) = 0.0d0 1585 !unita(2,3) = 0.0d0 1586 !unita(3,3) = 20.0d0 1587 end if 1588 value = beom_destroy(geom) 1589 if (.not. value) call errquit('cannot destroy geom',0, 1590 & GEOM_ERR) 1591 1592 end if 1593 1594 return 1595 end 1596 1597* *********************************** 1598* * * 1599* * control_find_box * 1600* * * 1601* *********************************** 1602 subroutine control_find_box(geom,box_type,orient,delta,unita) 1603 implicit none 1604 integer geom 1605 integer box_type 1606 logical orient 1607 real*8 delta 1608 real*8 unita(3,3) 1609 1610#include "bafdecls.fh" 1611#include "btdb.fh" 1612#include "beom.fh" 1613#include "errquit.fh" 1614 1615* *** local variables *** 1616 integer lwork 1617 parameter (lwork=9) 1618 double precision work(lwork) 1619 1620 integer ii,nion,ierr,i,j 1621 double precision q,rxyz(3),mtensor(3,3),frac(3),L1,L2,L3 1622 double precision x,y,z,meig(3) 1623 double precision Ixx,Iyx,Izx 1624 double precision Ixy,Iyy,Izy 1625 double precision Ixz,Iyz,Izz,Lmax,tmass 1626 double precision txy,txz,tyz,txx_yy,txx_zz,tyy_zz,tx,ty,tz 1627 1628 character*16 t 1629 logical value 1630 1631* **** external functions **** 1632 logical control_notqmmmq 1633 external control_notqmmmq 1634 1635 if (.not.geom_ncent(geom,nion)) then 1636 call errquit('cannot load nion from geom',0,GEOM_ERR) 1637 end if 1638 1639* *** find principle axes wrt origin (0,0,0) *** 1640 mtensor(1,1) = 1.0d0 1641 mtensor(2,1) = 0.0d0 1642 mtensor(3,1) = 0.0d0 1643 mtensor(1,2) = 0.0d0 1644 mtensor(2,2) = 1.0d0 1645 mtensor(3,2) = 0.0d0 1646 mtensor(1,3) = 0.0d0 1647 mtensor(2,3) = 0.0d0 1648 mtensor(3,3) = 1.0d0 1649 if (orient) then 1650 call dcopy(9,0.0d0,0,mtensor,1) 1651 txy = 0.0d0 1652 txz = 0.0d0 1653 tyz = 0.0d0 1654 txx_yy = 0.0d0 1655 txx_zz = 0.0d0 1656 tyy_zz = 0.0d0 1657 tx = 0.0d0 1658 ty = 0.0d0 1659 tz = 0.0d0 1660 do ii=1,nion 1661 value = geom_cent_get(geom,ii,t,rxyz,q) 1662 if (.not.control_notqmmmq(t)) then 1663 x = rxyz(1) 1664 y = rxyz(2) 1665 z = rxyz(3) 1666 1667 txy = txy + x*y 1668 txz = txz + x*z 1669 tyz = tyz + y*z 1670 1671 txx_yy = txx_yy + (x*x + y*y) 1672 txx_zz = txx_zz + (x*x + z*z) 1673 tyy_zz = tyy_zz + (y*y + z*z) 1674 1675 tx = tx + x 1676 ty = ty + y 1677 tz = tz + z 1678 end if 1679 end do 1680 tmass = dble(nion) 1681 Ixx = tyy_zz - ty*ty/tmass - tz*tz/tmass 1682 Iyy = txx_zz - tx*tx/tmass - tz*tz/tmass 1683 Izz = txx_yy - tx*tx/tmass - ty*ty/tmass 1684 Ixy = txy - tx*ty/tmass 1685 Iyx = Ixy 1686 Ixz = txz - tx*tz/tmass 1687 Izx = Ixz 1688 Iyz = tyz - ty*tz/tmass 1689 Izy = Iyz 1690 1691 mtensor(1,1) = Ixx 1692 mtensor(2,1) = -Ixy 1693 mtensor(3,1) = -Ixz 1694 1695 mtensor(1,2) = -Iyx 1696 mtensor(2,2) = Iyy 1697 mtensor(3,2) = -Iyz 1698 1699 mtensor(1,3) = -Izx 1700 mtensor(2,3) = -Izy 1701 mtensor(3,3) = Izz 1702 1703c **** longest dimension is along a1 **** 1704 ierr = 0 1705 call DSYEV('V','U',3,mtensor,3,meig,work,lwork,ierr) 1706 1707 1708 1709c !*** reorder eigenvectors - make longest dimesion along a3 **** 1710c x = mtensor(1,1) 1711c y = mtensor(2,1) 1712c z = mtensor(3,1) 1713c mtensor(1,1) = mtensor(1,3) 1714c mtensor(2,1) = mtensor(2,3) 1715c mtensor(3,1) = mtensor(3,3) 1716c mtensor(1,3) = x 1717c mtensor(2,3) = y 1718c mtensor(3,3) = z 1719 end if 1720 1721* ***************** 1722* *** cubic box *** 1723* ***************** 1724 if (box_type.eq.0) then 1725 1726* **** define L1 *** 1727 L1 = 0.0d0 1728 do ii=1,nion 1729 if (.not.geom_cent_get(geom,ii,t,rxyz,q)) then 1730 call errquit('cannot load center from geom',0,GEOM_ERR) 1731 end if 1732 if (.not.control_notqmmmq(t)) then 1733 frac(1) = mtensor(1,1)*rxyz(1) 1734 > + mtensor(2,1)*rxyz(2) 1735 > + mtensor(3,1)*rxyz(3) 1736 frac(2) = mtensor(1,2)*rxyz(1) 1737 > + mtensor(2,2)*rxyz(2) 1738 > + mtensor(3,2)*rxyz(3) 1739 frac(3) = mtensor(1,3)*rxyz(1) 1740 > + mtensor(2,3)*rxyz(2) 1741 > + mtensor(3,3)*rxyz(3) 1742 if ((frac(1)+delta).gt.( 0.5*L1)) L1 = 2.0d0*(frac(1)+delta) 1743 if ((frac(2)+delta).gt.( 0.5*L1)) L1 = 2.0d0*(frac(2)+delta) 1744 if ((frac(3)+delta).gt.( 0.5*L1)) L1 = 2.0d0*(frac(3)+delta) 1745 if ((frac(1)-delta).lt.(-0.5*L1)) 1746 > L1 = 2.0d0*(dabs(frac(1))+delta) 1747 if ((frac(2)-delta).lt.(-0.5*L1)) 1748 > L1 = 2.0d0*(dabs(frac(2))+delta) 1749 if ((frac(3)-delta).lt.(-0.5*L1)) 1750 > L1 = 2.0d0*(dabs(frac(3))+delta) 1751 end if 1752 end do 1753 !*** put threshold on smallest cubic box *** 1754 if (L1.lt.26.0d0) L1=26.0d0 1755 1756* **** define unit cell **** 1757 unita(1,1) = L1*mtensor(1,1) 1758 unita(2,1) = L1*mtensor(2,1) 1759 unita(3,1) = L1*mtensor(3,1) 1760 unita(1,2) = L1*mtensor(1,2) 1761 unita(2,2) = L1*mtensor(2,2) 1762 unita(3,2) = L1*mtensor(3,2) 1763 unita(1,3) = L1*mtensor(1,3) 1764 unita(2,3) = L1*mtensor(2,3) 1765 unita(3,3) = L1*mtensor(3,3) 1766 1767* ************************ 1768* *** orthorhombic box *** 1769* ************************ 1770 else if (box_type.eq.1) then 1771 1772* **** define L1, L2, and L3 **** 1773 L1 = 0.0d0 1774 L2 = 0.0d0 1775 L3 = 0.0d0 1776 do ii=1,nion 1777 if (.not.geom_cent_get(geom,ii,t,rxyz,q)) then 1778 call errquit('cannot load center from geom',0,GEOM_ERR) 1779 end if 1780 if (.not.control_notqmmmq(t)) then 1781 frac(1) = mtensor(1,1)*rxyz(1) 1782 > + mtensor(2,1)*rxyz(2) 1783 > + mtensor(3,1)*rxyz(3) 1784 frac(2) = mtensor(1,2)*rxyz(1) 1785 > + mtensor(2,2)*rxyz(2) 1786 > + mtensor(3,2)*rxyz(3) 1787 frac(3) = mtensor(1,3)*rxyz(1) 1788 > + mtensor(2,3)*rxyz(2) 1789 > + mtensor(3,3)*rxyz(3) 1790c write(*,*) "ii,frac=",ii,frac 1791c write(*,*) " =", 1792c > mtensor(1,1),mtensor(2,1),mtensor(3,1) 1793c write(*,*) " =",rxyz 1794c write(*,*) 1795 if ((frac(1)+delta).gt.(0.5d0*L1)) L1 = 2.0d0*(frac(1)+delta) 1796 if ((frac(2)+delta).gt.(0.5d0*L2)) L2 = 2.0d0*(frac(2)+delta) 1797 if ((frac(3)+delta).gt.(0.5d0*L3)) L3 = 2.0d0*(frac(3)+delta) 1798 if ((frac(1)-delta).lt.(-0.5d0*L1)) 1799 > L1 = 2.0d0*(dabs(frac(1))+delta) 1800 if ((frac(2)-delta).lt.(-0.5d0*L2)) 1801 > L2 = 2.0d0*(dabs(frac(2))+delta) 1802 if ((frac(3)-delta).lt.(-0.5d0*L3)) 1803 > L3 = 2.0d0*(dabs(frac(3))+delta) 1804 end if 1805 end do 1806 1807c write(*,*) "L1=",L1 1808c write(*,*) "L2=",L2 1809c write(*,*) "L3=",L3 1810 !*** put threshold on smallest box *** 1811 Lmax = L1 1812 if (Lmax.lt.L2) Lmax = L2 1813 if (Lmax.lt.L3) Lmax = L3 1814 if (Lmax.lt.26.0d0) then 1815 L1 = 19.0d0 1816 L2 = 19.0d0 1817 L3 = 19.0d0 1818 mtensor(1,1) = 1.0d0 1819 mtensor(2,1) = 1.0d0 1820 mtensor(3,1) = 0.0d0 1821 mtensor(1,2) = 1.0d0 1822 mtensor(2,2) = 0.0d0 1823 mtensor(3,2) = 1.0d0 1824 mtensor(1,3) = 0.0d0 1825 mtensor(2,3) = 1.0d0 1826 mtensor(3,3) = 1.0d0 1827 end if 1828 1829* **** define unit cell **** 1830 unita(1,1) = L1*mtensor(1,1) 1831 unita(2,1) = L1*mtensor(2,1) 1832 unita(3,1) = L1*mtensor(3,1) 1833 unita(1,2) = L2*mtensor(1,2) 1834 unita(2,2) = L2*mtensor(2,2) 1835 unita(3,2) = L2*mtensor(3,2) 1836 unita(1,3) = L3*mtensor(1,3) 1837 unita(2,3) = L3*mtensor(2,3) 1838 unita(3,3) = L3*mtensor(3,3) 1839 1840* **** unknown box type **** 1841 else 1842 call errquit('invalid box_type',0,0) 1843 end if 1844 1845 1846 return 1847 end 1848 1849* *************************** 1850* * * 1851* * control_notqmmmq * 1852* * * 1853* *************************** 1854 logical function control_notqmmmq(string) 1855 implicit none 1856 character*16 string 1857 1858 logical qmmmq 1859 1860 qmmmq = .false. 1861 if (index(string,'^').gt.0) qmmmq = .true. 1862 if (index(string,'x').eq.1) qmmmq = .true. 1863 if (index(string,'X').eq.1) qmmmq = .true. 1864 if (index(string,'bq').eq.1) qmmmq = .true. 1865 if (index(string,'Bq').eq.1) qmmmq = .true. 1866 if (index(string,'bQ').eq.1) qmmmq = .true. 1867 if (index(string,'BQ').eq.1) qmmmq = .true. 1868 1869 1870 control_notqmmmq = qmmmq 1871 return 1872 end 1873 1874 1875 1876* *********************************** 1877* * * 1878* * control_move * 1879* * * 1880* *********************************** 1881 logical function control_move() 1882 implicit none 1883 1884#include "control.fh" 1885 1886 control_move = move 1887 return 1888 end 1889 1890 1891* *********************************** 1892* * * 1893* * control_rotation * 1894* * * 1895* *********************************** 1896 logical function control_rotation() 1897 implicit none 1898 1899#include "control.fh" 1900 1901 control_rotation = rotation 1902 return 1903 end 1904 1905* *********************************** 1906* * * 1907* * control_dof_rotation * 1908* * * 1909* *********************************** 1910 logical function control_dof_rotation() 1911 implicit none 1912 1913#include "control.fh" 1914 1915 control_dof_rotation = dof_rotation 1916 return 1917 end 1918 1919* *********************************** 1920* * * 1921* * control_translation * 1922* * * 1923* *********************************** 1924 logical function control_translation() 1925 implicit none 1926 1927#include "control.fh" 1928 1929 control_translation = translation 1930 return 1931 end 1932 1933* *********************************** 1934* * * 1935* * control_dof_translation * 1936* * * 1937* *********************************** 1938 logical function control_dof_translation() 1939 implicit none 1940 1941#include "control.fh" 1942 1943 control_dof_translation = dof_translation 1944 return 1945 end 1946 1947 1948 1949 1950 1951* *********************************** 1952* * * 1953* * control_translate_vector * 1954* * * 1955* *********************************** 1956 subroutine control_translate_vector(rtrans) 1957 implicit none 1958 real*8 rtrans(3) 1959 1960#include "bafdecls.fh" 1961#include "btdb.fh" 1962 1963* **** control_rtdb common block **** 1964 integer rtdb 1965 common / control_rtdb1 / rtdb 1966 1967 if (.not.btdb_get(rtdb,'nwpw:translate_vector', 1968 > mt_dbl,3,rtrans)) then 1969 rtrans(1) = 0.0d0 1970 rtrans(2) = 0.0d0 1971 rtrans(3) = 0.0d0 1972 end if 1973 return 1974 end 1975 1976* *********************************** 1977* * * 1978* * control_translate_reorder * 1979* * * 1980* *********************************** 1981 subroutine control_translate_reorder(reorder) 1982 implicit none 1983 logical reorder 1984 1985#include "bafdecls.fh" 1986#include "btdb.fh" 1987 1988* **** control_rtdb common block **** 1989 integer rtdb 1990 common / control_rtdb1 / rtdb 1991 1992 if (.not.btdb_get(rtdb,'nwpw:translate_reorder', 1993 > mt_log,1,reorder)) then 1994 reorder = .true. 1995 end if 1996 return 1997 end 1998 1999 2000 2001* *********************************** 2002* * * 2003* * control_translate_geom_name * 2004* * * 2005* *********************************** 2006 subroutine control_translate_geom_name(geom_name) 2007 implicit none 2008 character*(*) geom_name 2009 2010#include "bafdecls.fh" 2011#include "btdb.fh" 2012 2013* **** control_rtdb common block **** 2014 integer rtdb 2015 common / control_rtdb1 / rtdb 2016 2017 2018 if (.not.btdb_cget(rtdb,'nwpw:translate_geom_name', 2019 > 1,geom_name)) then 2020 geom_name = "translated_geometry" 2021 end if 2022 return 2023 end 2024 2025 2026 2027 2028* *********************************** 2029* * * 2030* * control_out_of_time * 2031* * * 2032* *********************************** 2033 2034* This function is used to estimate if there is 2035* enough time to perform another iteration. The 2036* routine control_read intializes this routine. To 2037* determine if there is enough time left to do another 2038* iteration this routine uses estimates for the amount 2039* of time to finish the simulation (est_finish_time) and 2040* the amount of time to perform another iteration step 2041* (est_step_time). Where 2042* 2043* est_finish_time = 2*(time elapsed from call to control_read 2044* to the first call to control_out_of_time) 2045* 2046* est_step_time = (time elapsed between successive calls to 2047* control_out_of_time) 2048* 2049* Uses: control_blktime common block located in control.fh 2050* util_test_time_remaining 2051* created: 5-8-2002 2052 2053 logical function control_out_of_time() 2054 implicit none 2055 2056#include "control.fh" 2057 2058* **** control_rtdb common block **** 2059 integer rtdb 2060 common / control_rtdb1 / rtdb 2061 2062* **** local variables **** 2063 logical value 2064 integer required_time 2065 2066* **** external functions **** 2067 logical util_test_time_remaining 2068 external util_test_time_remaining 2069 2070 call current_second(cpu2_time) 2071 2072* **** This is the first time this routine has been called **** 2073 if (est_finish_time.eq.-1) then 2074 est_finish_time = 2*int(cpu2_time-cpu1_time) ! crude estimate 2075 value = .false. 2076 2077* **** This routine has been called two or more times **** 2078 else 2079 est_step_time = int(cpu2_time-cpu1_time)+1 ! no statistical info used 2080 required_time = est_step_time + est_finish_time 2081 value = .not.util_test_time_remaining(rtdb,required_time) 2082 end if 2083 2084 cpu1_time = cpu2_time 2085 2086 control_out_of_time = value 2087 return 2088 end 2089 2090 2091* *********************************** 2092* * * 2093* * control_frac_coord * 2094* * * 2095* *********************************** 2096 logical function control_frac_coord() 2097 implicit none 2098 2099#include "control.fh" 2100 2101 control_frac_coord = frac_coord 2102 return 2103 end 2104 2105 2106 2107 2108 2109* *********************************** 2110* * * 2111* * control_code * 2112* * * 2113* *********************************** 2114 integer function control_code() 2115 implicit none 2116 2117#include "control.fh" 2118 2119 control_code = code 2120 return 2121 end 2122 2123 2124 2125* *********************************** 2126* * * 2127* * control_ngrid * 2128* * * 2129* *********************************** 2130 integer function control_ngrid(ijk) 2131 implicit none 2132 integer ijk 2133 2134#include "control.fh" 2135 2136 control_ngrid = ngrid(ijk) 2137 return 2138 end 2139 2140 2141* *********************************** 2142* * * 2143* * control_ngrid_small * 2144* * * 2145* *********************************** 2146 integer function control_ngrid_small(ijk) 2147 implicit none 2148 integer ijk 2149 2150#include "control.fh" 2151 2152 control_ngrid_small = ngrid_small(ijk) 2153 return 2154 end 2155 2156* *********************************** 2157* * * 2158* * control_has_ngrid_small * 2159* * * 2160* *********************************** 2161 logical function control_has_ngrid_small() 2162 implicit none 2163 2164#include "control.fh" 2165 2166 control_has_ngrid_small = has_ngrid_small 2167 return 2168 end 2169 2170 2171* *********************************** 2172* * * 2173* * control_it_in * 2174* * * 2175* *********************************** 2176 integer function control_it_in() 2177 implicit none 2178 2179#include "control.fh" 2180 2181 control_it_in = loop(1) 2182 return 2183 end 2184 2185 2186* *********************************** 2187* * * 2188* * control_it_out * 2189* * * 2190* *********************************** 2191 integer function control_it_out() 2192 implicit none 2193 2194#include "control.fh" 2195 2196 control_it_out = loop(2) 2197 return 2198 end 2199 2200* *********************************** 2201* * * 2202* * control_bo_steps_in * 2203* * * 2204* *********************************** 2205 integer function control_bo_steps_in() 2206 implicit none 2207 2208#include "control.fh" 2209 2210 control_bo_steps_in = bo_steps(1) 2211 return 2212 end 2213 2214* *********************************** 2215* * * 2216* * control_bo_steps_out * 2217* * * 2218* *********************************** 2219 integer function control_bo_steps_out() 2220 implicit none 2221 2222#include "control.fh" 2223 2224 control_bo_steps_out = bo_steps(2) 2225 return 2226 end 2227 2228* *********************************** 2229* * * 2230* * control_bo_algorithm * 2231* * * 2232* *********************************** 2233 integer function control_bo_algorithm() 2234 implicit none 2235 2236#include "control.fh" 2237 2238 control_bo_algorithm = bo_algorithm 2239 return 2240 end 2241 2242 2243 2244* *********************************** 2245* * * 2246* * control_time_step * 2247* * * 2248* *********************************** 2249 real*8 function control_time_step() 2250 implicit none 2251 2252#include "control.fh" 2253 2254 control_time_step = time_step 2255 return 2256 end 2257 2258* *********************************** 2259* * * 2260* * control_bo_time_step * 2261* * * 2262* *********************************** 2263 real*8 function control_bo_time_step() 2264 implicit none 2265 2266#include "control.fh" 2267 2268 control_bo_time_step = bo_time_step 2269 2270 return 2271 end 2272 2273* *********************************** 2274* * * 2275* * control_ion_time_step * 2276* * * 2277* *********************************** 2278* Used by ion.F, the ion_time_step is 2279* set to time_step if Car-Parrinello 2280* and it is set to bo_time_step if 2281* Born-Oppenheimer. 2282 real*8 function control_ion_time_step() 2283 implicit none 2284 2285#include "control.fh" 2286 2287 !*** BO dynamics **** 2288 if (code.eq.11) then 2289 control_ion_time_step = bo_time_step 2290 2291 !*** CP dynamics **** 2292 else 2293 control_ion_time_step = time_step 2294 end if 2295 2296 return 2297 end 2298 2299 2300 2301* *********************************** 2302* * * 2303* * control_fake_mass * 2304* * * 2305* *********************************** 2306 real*8 function control_fake_mass() 2307 implicit none 2308 2309#include "control.fh" 2310 2311 control_fake_mass = fake_mass 2312 return 2313 end 2314 2315 2316* *********************************** 2317* * * 2318* * control_bo_fake_mass * 2319* * * 2320* *********************************** 2321 real*8 function control_bo_fake_mass() 2322 implicit none 2323 2324#include "control.fh" 2325 2326 control_bo_fake_mass = bo_fake_mass 2327 return 2328 end 2329 2330* *********************************** 2331* * * 2332* * control_ks_alpha * 2333* * * 2334* *********************************** 2335 real*8 function control_ks_alpha() 2336 implicit none 2337 2338#include "control.fh" 2339 2340 control_ks_alpha = ks_alpha 2341 return 2342 end 2343 2344* *********************************** 2345* * * 2346* * control_fractional_alpha * 2347* * * 2348* *********************************** 2349 real*8 function control_fractional_alpha() 2350 implicit none 2351 2352#include "control.fh" 2353 2354 control_fractional_alpha = fractional_alpha 2355 return 2356 end 2357 2358 2359* *********************************** 2360* * * 2361* * control_kerker_g0 * 2362* * * 2363* *********************************** 2364 real*8 function control_kerker_g0() 2365 implicit none 2366 2367#include "control.fh" 2368 2369 control_kerker_g0 = kerker_g0 2370 return 2371 end 2372 2373 2374* *********************************** 2375* * * 2376* * control_scf_algorithm * 2377* * * 2378* *********************************** 2379 integer function control_scf_algorithm() 2380 implicit none 2381 2382#include "control.fh" 2383 2384 control_scf_algorithm = scf_algorithm 2385 return 2386 end 2387 2388* *********************************** 2389* * * 2390* * control_diis_histories * 2391* * * 2392* *********************************** 2393 integer function control_diis_histories() 2394 implicit none 2395 2396#include "control.fh" 2397 2398 control_diis_histories = diis_histories 2399 return 2400 end 2401 2402* *********************************** 2403* * * 2404* * control_ks_algorithm * 2405* * * 2406* *********************************** 2407 integer function control_ks_algorithm() 2408 implicit none 2409 2410#include "control.fh" 2411 2412 control_ks_algorithm = ks_algorithm 2413 return 2414 end 2415 2416 2417* *********************************** 2418* * * 2419* * control_ks_maxit_orb * 2420* * * 2421* *********************************** 2422 integer function control_ks_maxit_orb() 2423 implicit none 2424 2425#include "control.fh" 2426 2427 control_ks_maxit_orb = maxit_orb 2428 return 2429 end 2430 2431* *********************************** 2432* * * 2433* * control_ks_maxit_orbs * 2434* * * 2435* *********************************** 2436 integer function control_ks_maxit_orbs() 2437 implicit none 2438 2439#include "control.fh" 2440 2441 control_ks_maxit_orbs = maxit_orbs 2442 return 2443 end 2444 2445 2446 2447* *********************************** 2448* * * 2449* * control_tole * 2450* * * 2451* *********************************** 2452 real*8 function control_tole() 2453 implicit none 2454 2455#include "control.fh" 2456 2457 control_tole = tolerances(1) 2458 return 2459 end 2460 2461 2462* *********************************** 2463* * * 2464* * control_tolc * 2465* * * 2466* *********************************** 2467 real*8 function control_tolc() 2468 implicit none 2469 2470#include "control.fh" 2471 2472 control_tolc = tolerances(2) 2473 return 2474 end 2475 2476 2477* *********************************** 2478* * * 2479* * control_tolr * 2480* * * 2481* *********************************** 2482 real*8 function control_tolr() 2483 implicit none 2484 2485#include "control.fh" 2486 2487 control_tolr = tolerances(3) 2488 return 2489 end 2490 2491* *********************************** 2492* * * 2493* * control_rte * 2494* * * 2495* *********************************** 2496 real*8 function control_rte() 2497 implicit none 2498 2499#include "control.fh" 2500 2501 control_rte = scaling(1) 2502 return 2503 end 2504 2505* *********************************** 2506* * * 2507* * control_rti * 2508* * * 2509* *********************************** 2510 real*8 function control_rti() 2511 implicit none 2512 2513#include "control.fh" 2514 2515 control_rti = scaling(2) 2516 return 2517 end 2518 2519 2520* *********************************** 2521* * * 2522* * control_unita * 2523* * * 2524* *********************************** 2525 real*8 function control_unita(i,j) 2526 implicit none 2527 integer i,j 2528 2529#include "control.fh" 2530 2531 control_unita = unita(i,j) 2532 return 2533 end 2534 2535 2536* *********************************** 2537* * * 2538* * control_set_unita * 2539* * * 2540* *********************************** 2541 subroutine control_set_unita(unita_in) 2542 implicit none 2543 real*8 unita_in(*) 2544 2545#include "bafdecls.fh" 2546#include "btdb.fh" 2547#include "errquit.fh" 2548#include "control.fh" 2549 2550* **** local variables **** 2551 integer taskid,MASTER 2552 parameter (MASTER=0) 2553 logical value 2554 integer l 2555 character*50 rtdb_unita 2556 2557* **** control_rtdb common block **** 2558 integer rtdb 2559 common / control_rtdb1 / rtdb 2560 2561 2562 call dcopy(9,unita_in,1,unita,1) 2563 2564* **** put unita on the rtdb **** 2565 l = index(cell_name,' ') - 1 2566 rtdb_unita = cell_name(1:l)//':unita' 2567 2568 call Parallel_taskid(taskid) 2569 value=btdb_parallel(.false.) 2570 if (taskid.eq.MASTER) then 2571 if (.not.btdb_put(rtdb,rtdb_unita,mt_dbl,9,unita)) 2572 > call errquit('control_set_unita:writing unita',0,RTDB_ERR) 2573 end if 2574 value=btdb_parallel(.true.) 2575 return 2576 end 2577 2578 2579* *********************************** 2580* * * 2581* * control_get_unita * 2582* * * 2583* *********************************** 2584 subroutine control_get_unita(unita_out) 2585 implicit none 2586 real*8 unita_out(*) 2587 2588#include "control.fh" 2589 2590 call dcopy(9,unita,1,unita_out,1) 2591 return 2592 end 2593 2594 2595 2596* *********************************** 2597* * * 2598* * control_set_unita_frozen * 2599* * * 2600* *********************************** 2601 subroutine control_set_unita_frozen(unita_in) 2602 implicit none 2603 real*8 unita_in(*) 2604 2605#include "bafdecls.fh" 2606#include "btdb.fh" 2607#include "errquit.fh" 2608#include "control.fh" 2609 2610* **** local variables **** 2611 integer taskid,MASTER 2612 parameter (MASTER=0) 2613 logical value 2614 integer l 2615 character*50 rtdb_unita 2616 2617* **** control_rtdb common block **** 2618 integer rtdb 2619 common / control_rtdb1 / rtdb 2620 2621 call dcopy(9,unita_in,1,unita_frozen,1) 2622 2623* **** put unita_frozen on the rtdb **** 2624 l = index(cell_name,' ') - 1 2625 rtdb_unita = cell_name(1:l)//':unita_frozen' 2626 2627 call Parallel_taskid(taskid) 2628 value=btdb_parallel(.false.) 2629 if (taskid.eq.MASTER) then 2630 if (.not.btdb_put(rtdb,rtdb_unita,mt_dbl,9,unita)) 2631 > call errquit('control_set_unita:writing unita',0,RTDB_ERR) 2632 end if 2633 value=btdb_parallel(.true.) 2634 return 2635 end 2636 2637 2638* *********************************** 2639* * * 2640* * control_unita_frozen * 2641* * * 2642* *********************************** 2643 real*8 function control_unita_frozen(i,j) 2644 implicit none 2645 integer i,j 2646 2647#include "control.fh" 2648 2649 control_unita_frozen = unita_frozen(i,j) 2650 return 2651 end 2652 2653* *********************************** 2654* * * 2655* * control_frozen * 2656* * * 2657* *********************************** 2658 logical function control_frozen() 2659 implicit none 2660 2661#include "control.fh" 2662 2663 control_frozen = frozen 2664 return 2665 end 2666 2667 2668* *********************************** 2669* * * 2670* * control_boundry * 2671* * * 2672* *********************************** 2673 character*12 function control_boundry() 2674 implicit none 2675 2676#include "control.fh" 2677 2678 control_boundry = boundry 2679 return 2680 end 2681 2682 2683c* *********************************** 2684c* * * 2685c* * control_pspnames * 2686c* * * 2687c* *********************************** 2688c character*20 function control_pspnames(i) 2689c implicit none 2690c integer i 2691c 2692c#include "control.fh" 2693c 2694c control_pspnames = pspnames(i) 2695c return 2696c end 2697c 2698c 2699c* *********************************** 2700c* * * 2701c* * control_pspstressnames * 2702c* * * 2703c* *********************************** 2704c character*20 function control_pspstressnames(i) 2705c implicit none 2706c integer i 2707c 2708c integer ind 2709c character*20 pspname 2710c character*20 control_pspnames 2711c external control_pspnames 2712c 2713c pspname = control_pspnames(i) 2714c ind = index(pspname,' ') -1 2715c pspname = pspname(1:ind)//'2' 2716c 2717c control_pspstressnames = pspname 2718c return 2719c end 2720 2721* *********************************** 2722* * * 2723* * control_npsp * 2724* * * 2725* *********************************** 2726 integer function control_npsp() 2727 implicit none 2728 2729#include "control.fh" 2730 2731 control_npsp = npsp 2732 return 2733 end 2734 2735 2736 2737* *********************************** 2738* * * 2739* * control_ecut * 2740* * * 2741* *********************************** 2742 real*8 function control_ecut() 2743 implicit none 2744 2745#include "control.fh" 2746 control_ecut = ecut 2747 return 2748 end 2749 2750 2751 2752* *********************************** 2753* * * 2754* * control_wcut * 2755* * * 2756* *********************************** 2757 real*8 function control_wcut() 2758 implicit none 2759 2760#include "control.fh" 2761 2762 control_wcut = wcut 2763 return 2764 end 2765 2766 2767* *********************************** 2768* * * 2769* * control_rcut * 2770* * * 2771* *********************************** 2772 real*8 function control_rcut() 2773 implicit none 2774 2775#include "control.fh" 2776 2777 control_rcut = rcut 2778 return 2779 end 2780 2781* *********************************** 2782* * * 2783* * control_ncut * 2784* * * 2785* *********************************** 2786 integer function control_ncut() 2787 implicit none 2788 2789#include "control.fh" 2790 2791 control_ncut = ncut 2792 return 2793 end 2794 2795 2796 2797 2798* *********************************** 2799* * * 2800* * control_output_psi * 2801* * * 2802* *********************************** 2803 character*50 function control_output_psi() 2804 implicit none 2805 2806#include "control.fh" 2807 2808 control_output_psi = output_wavefunction_filename 2809 return 2810 end 2811 2812* *********************************** 2813* * * 2814* * control_output_epsi * 2815* * * 2816* *********************************** 2817 character*50 function control_output_epsi() 2818 implicit none 2819 2820#include "control.fh" 2821 2822 control_output_epsi = output_ewavefunction_filename 2823 return 2824 end 2825 2826 2827* *********************************** 2828* * * 2829* * control_input_psi * 2830* * * 2831* *********************************** 2832 character*50 function control_input_psi() 2833 implicit none 2834 2835#include "control.fh" 2836 2837 control_input_psi = input_wavefunction_filename 2838 return 2839 end 2840 2841 2842* *********************************** 2843* * * 2844* * control_input_epsi * 2845* * * 2846* *********************************** 2847 character*50 function control_input_epsi() 2848 implicit none 2849 2850#include "control.fh" 2851 2852 control_input_epsi = input_ewavefunction_filename 2853 return 2854 end 2855 2856 2857* *********************************** 2858* * * 2859* * control_output_v_psi * 2860* * * 2861* *********************************** 2862 character*50 function control_output_v_psi() 2863 implicit none 2864 2865#include "control.fh" 2866 2867 control_output_v_psi = output_v_wavefunction_filename 2868 return 2869 end 2870 2871 2872* *********************************** 2873* * * 2874* * control_input_v_psi * 2875* * * 2876* *********************************** 2877 character*50 function control_input_v_psi() 2878 implicit none 2879 2880#include "control.fh" 2881 2882 control_input_v_psi = input_v_wavefunction_filename 2883 return 2884 end 2885 2886 2887* *********************************** 2888* * * 2889* * control_use_fractional_rho * 2890* * * 2891* *********************************** 2892 logical function control_use_fractional_rho() 2893 implicit none 2894 2895#include "bafdecls.fh" 2896#include "btdb.fh" 2897#include "errquit.fh" 2898#include "control.fh" 2899 2900* **** control_rtdb common block **** 2901 integer rtdb 2902 common / control_rtdb1 / rtdb 2903 2904 logical use_rho 2905 2906 if (.not.btdb_get(rtdb,"nwpw:use_fractional_rho", 2907 > mt_log,1,use_rho)) then 2908 use_rho = .true. 2909 end if 2910 control_use_fractional_rho = use_rho 2911 return 2912 end 2913 2914 2915* *********************************** 2916* * * 2917* * control_input_rho * 2918* * * 2919* *********************************** 2920 2921 character*50 function control_input_rho() 2922 implicit none 2923 2924#include "bafdecls.fh" 2925#include "btdb.fh" 2926#include "errquit.fh" 2927#include "control.fh" 2928 2929* **** control_rtdb common block **** 2930 integer rtdb 2931 common / control_rtdb1 / rtdb 2932 2933 character*50 filename 2934 2935 if (.not.btdb_cget(rtdb,'nwpw:input rho',1,filename)) then 2936 call util_file_prefix('rho',filename) 2937 end if 2938 control_input_rho = filename 2939 return 2940 end 2941 2942* *********************************** 2943* * * 2944* * control_output_rho * 2945* * * 2946* *********************************** 2947 2948 character*50 function control_output_rho() 2949 implicit none 2950 2951#include "bafdecls.fh" 2952#include "btdb.fh" 2953#include "errquit.fh" 2954#include "control.fh" 2955 2956* **** control_rtdb common block **** 2957 integer rtdb 2958 common / control_rtdb1 / rtdb 2959 2960 character*50 filename 2961 2962 if (.not.btdb_cget(rtdb,'nwpw:output rho',1,filename)) then 2963 call util_file_prefix('rho',filename) 2964 end if 2965 control_output_rho = filename 2966 return 2967 end 2968 2969 2970 2971 2972* *********************************** 2973* * * 2974* * control_xyz * 2975* * * 2976* *********************************** 2977 character*50 function control_xyz() 2978 implicit none 2979 2980 2981#include "control.fh" 2982 2983 control_xyz = xyz_filename 2984 return 2985 end 2986 2987* *********************************** 2988* * * 2989* * control_cell_name * 2990* * * 2991* *********************************** 2992 character*50 function control_cell_name() 2993 implicit none 2994 2995 2996#include "control.fh" 2997 2998 control_cell_name = cell_name 2999 return 3000 end 3001 3002 3003 3004* *********************************** 3005* * * 3006* * control_gga * 3007* * * 3008* *********************************** 3009 integer function control_gga() 3010 implicit none 3011 3012#include "control.fh" 3013 3014 control_gga = gga 3015 return 3016 end 3017 3018* *********************************** 3019* * * 3020* * control_only_lda * 3021* * * 3022* *********************************** 3023 logical function control_only_lda() 3024 implicit none 3025 3026#include "control.fh" 3027 3028* **** external functions **** 3029 logical nwpwxc_is_on,nwpwxc_is_lda 3030 external nwpwxc_is_on,nwpwxc_is_lda 3031 3032 control_only_lda = (.not.nwpwxc_is_on().and.gga.eq.0).or. 3033 > (nwpwxc_is_on().and.nwpwxc_is_lda()) 3034 return 3035 end 3036 3037 3038 3039* *********************************** 3040* * * 3041* * control_is_grimme2 * 3042* * * 3043* *********************************** 3044 logical function control_is_grimme2() 3045 implicit none 3046 3047#include "control.fh" 3048 3049 control_is_grimme2 = is_grimme2 3050 return 3051 end 3052 3053 3054* *********************************** 3055* * * 3056* * control_has_disp * 3057* * * 3058* *********************************** 3059 logical function control_has_disp() 3060 implicit none 3061 3062#include "control.fh" 3063 3064 control_has_disp = has_disp 3065 return 3066 end 3067 3068 3069* *********************************** 3070* * * 3071* * control_options_disp * 3072* * * 3073* *********************************** 3074 character*80 function control_options_disp() 3075 implicit none 3076 3077#include "control.fh" 3078 3079 control_options_disp = options_disp 3080 return 3081 end 3082 3083* *********************************** 3084* * * 3085* * control_has_vdw * 3086* * * 3087* *********************************** 3088 logical function control_has_vdw() 3089 implicit none 3090 3091#include "control.fh" 3092 3093 control_has_vdw = has_vdw 3094 return 3095 end 3096 3097* *********************************** 3098* * * 3099* * control_is_vdw2 * 3100* * * 3101* *********************************** 3102 logical function control_is_vdw2() 3103 implicit none 3104 3105#include "control.fh" 3106 3107 control_is_vdw2 = is_vdw2 3108 return 3109 end 3110 3111 3112* *********************************** 3113* * * 3114* * control_multiplicity * 3115* * * 3116* *********************************** 3117 integer function control_multiplicity() 3118 implicit none 3119 3120#include "control.fh" 3121 3122 control_multiplicity = multiplicity 3123 return 3124 end 3125 3126* *********************************** 3127* * * 3128* * control_multiplicity_set * 3129* * * 3130* *********************************** 3131 subroutine control_multiplicity_set(new_multiplicity) 3132 implicit none 3133 integer new_multiplicity 3134 3135#include "control.fh" 3136 3137 multiplicity = new_multiplicity 3138 return 3139 end 3140 3141* ***************************************** 3142* * * 3143* * control_check_charge_multiplicity * 3144* * * 3145* ***************************************** 3146 logical function control_check_charge_multiplicity() 3147 implicit none 3148 3149#include "control.fh" 3150 3151* *** local variables *** 3152 logical check 3153 real*8 icharge,tcharge,t 3154 integer mult,x,x_wf,nextra_orbs 3155 integer ispin_wf,ne_wf(2),x_f 3156 3157* ***** local functions **** 3158 logical psi_filefind 3159 external psi_filefind 3160 real*8 control_TotalCharge 3161 external control_TotalCharge 3162 real*8 ion_TotalCharge_qm 3163 external ion_TotalCharge_qm 3164 !integer control_frac_occ_extra_orbitals 3165 !external control_frac_occ_extra_orbitals 3166 integer control_fractional_orbitals 3167 external control_fractional_orbitals 3168 logical control_mult_fixed 3169 external control_mult_fixed 3170 3171 3172 3173* **** check wavefunction file **** 3174 if (psi_filefind()) then 3175 3176* **** get mult and e-charge from wavefunction file **** 3177 call psi_get_ne(ispin_wf,ne_wf) 3178 3179c write(*,*) "ne_wf = ",ne_wf 3180c x_wf0 = ne_wf(1)+ne_wf(2) 3181c mult0 = ne_wf(1)-ne_wf(2) + 1 3182 !nextra_orbs = control_frac_occ_extra_orbitals() 3183 if (ispin_wf.eq.1) then 3184 ne_wf(1) = ne_wf(1) - control_fractional_orbitals(1) 3185 else 3186 ne_wf(1) = ne_wf(1) - control_fractional_orbitals(1) 3187 ne_wf(2) = ne_wf(2) - control_fractional_orbitals(2) 3188 end if 3189c write(*,*) "after ne_wf = ",ne_wf 3190 3191 x_wf = ne_wf(1)+ne_wf(2) 3192 mult = ne_wf(1)-ne_wf(2) + 1 3193 if (ispin_wf.eq.1) then 3194c x_wf0 = 2*x_wf0 3195c mult0 = 1 3196 x_wf = 2*x_wf 3197 mult = 1 3198 end if 3199 3200* **** get mult and e-charge from control **** 3201 tcharge = control_TotalCharge() 3202 icharge = ion_TotalCharge_qm() 3203 t = icharge - tcharge !** total number of electrons ** 3204 x = int(NINT(t)) 3205 3206 if (.not.control_mult_fixed()) then 3207 3208* **** reassign spin to agree with total number of electrons **** 3209 if ((mod(x,2).ne.0).and.(ispin.eq.1)) then !** odd number of electrons ** 3210 ispin = 2 3211 end if 3212 3213* **** reassign multiplicity to agree with total number of electrons **** 3214* *** odd number of electrons and mult odd *** 3215 if ((mod(x,2).ne.0) .and.(mod(multiplicity,2).ne.0)) then 3216 multiplicity = multiplicity - 1 3217 do while (multiplicity.gt.(x+1)) 3218 multiplicity = multiplicity - 2 3219 end do 3220 if (multiplicity.lt.1) multiplicity = 2 3221 end if 3222* *** even number of electrons and mult even *** 3223 if ((mod(x,2).eq.0) .and.(mod(multiplicity,2).eq.0)) then 3224 multiplicity = multiplicity - 1 3225 do while (multiplicity.gt.(x+1)) 3226 multiplicity = multiplicity - 2 3227 end do 3228 if (multiplicity.lt.1) multiplicity = 1 3229 end if 3230 3231 3232* **** compare multiplicity, charge, and ispin **** 3233c check = ( (((mult.eq.multiplicity).and.(x_wf.eq.x)) .or. 3234c > ((mult0.eq.multiplicity).and.(x_wf0.eq.x))).and. 3235c > (ispin.eq.ispin_wf)) 3236 check = ( (mult.eq.multiplicity).and.(x_wf.eq.x).and. 3237 > (ispin.eq.ispin_wf)) 3238 3239 if (ispin_wf.eq.3) then 3240 check=((ne_wf(1).eq.ne_wf(2)).and. 3241 > (x.eq.ne_wf(1) ) ) 3242 end if 3243 3244 else 3245* **** compare multiplicity, and ispin **** 3246 check = ((mult.eq.multiplicity).and.(ispin.eq.ispin_wf)) 3247 end if 3248* **** no wavefunction file *** 3249 else 3250 check = .false. 3251 end if 3252 3253 control_check_charge_multiplicity = check 3254 return 3255 end 3256 3257 3258 3259* ***************************************** 3260* * * 3261* * control_check_number_virtuals * 3262* * * 3263* ***************************************** 3264 logical function control_check_number_virtuals() 3265 implicit none 3266 3267#include "control.fh" 3268 3269* *** local variables *** 3270 logical check 3271 integer ispin_wf,ne_wf(2),ne(2) 3272 3273* ***** local functions **** 3274 logical epsi_filefind 3275 external epsi_filefind 3276 integer control_excited_ne 3277 external control_excited_ne 3278 3279* **** check wavefunction file **** 3280 if (epsi_filefind()) then 3281 3282* **** get mult and e-charge from wavefunction file **** 3283 call psi_get_ne_excited(ispin_wf,ne_wf) 3284 ne(1) = 0 3285 ne(2) = 0 3286 ne(1) = control_excited_ne(1) 3287 if (ispin.eq.2) ne(2) = control_excited_ne(2) 3288 3289 check = ((ne(1).eq.ne_wf(1)).and. 3290 > (ne(2).eq.ne_wf(2)).and. 3291 > (ispin.eq.ispin_wf)) 3292 3293* **** no wavefunction file *** 3294 else 3295 check = .false. 3296 end if 3297 3298 control_check_number_virtuals = check 3299 return 3300 end 3301 3302 3303 3304* *********************************** 3305* * * 3306* * control_ispin * 3307* * * 3308* *********************************** 3309 integer function control_ispin() 3310 implicit none 3311 3312#include "control.fh" 3313 3314 control_ispin = ispin 3315 return 3316 end 3317 3318* *********************************** 3319* * * 3320* * control_ispin_set * 3321* * * 3322* *********************************** 3323 subroutine control_ispin_set(new_ispin) 3324 implicit none 3325 integer new_ispin 3326 3327#include "control.fh" 3328 3329 ispin = new_ispin 3330 return 3331 end 3332 3333 3334* ******************************************* 3335* * * 3336* * control_gradient_iterations * 3337* * * 3338* ******************************************* 3339 subroutine control_gradient_iterations() 3340 implicit none 3341 3342#include "control.fh" 3343 3344 loop(1) = 1 3345 loop(2) = 1 3346 3347 return 3348 end 3349 3350* *********************************** 3351* * * 3352* * control_version * 3353* * * 3354* *********************************** 3355 integer function control_version() 3356 implicit none 3357 3358#include "inp.fh" 3359#include "control.fh" 3360 3361* **** local variables **** 3362 integer l,version 3363 3364 l =index(boundry,' ') - 1 3365 3366 version = 3 3367 if (inp_compare(.false.,boundry(1:l),'periodic')) version=3 3368 if (inp_compare(.false.,boundry(1:l),'aperiodic')) version=4 3369 3370 control_version = version 3371 return 3372 end 3373 3374 3375* ************************ 3376* * * 3377* * control_Nose * 3378* * * 3379* ************************ 3380 logical function control_Nose() 3381 implicit none 3382 3383* **** control_nose common block **** 3384 logical nose 3385 real*8 Pe,Te,Pr,Tr 3386 common / control_nblock / Pe,Te,Pr,Tr,nose 3387 3388 3389 control_Nose = nose 3390 return 3391 end 3392 3393 3394* **************************** 3395* * * 3396* * control_Nose_Pe * 3397* * * 3398* **************************** 3399 real*8 function control_Nose_Pe() 3400 implicit none 3401 3402* **** control_nose common block **** 3403 logical nose 3404 real*8 Pe,Te,Pr,Tr 3405 common / control_nblock / Pe,Te,Pr,Tr,nose 3406 3407 3408 control_Nose_Pe = Pe 3409 return 3410 end 3411 3412 3413* **************************** 3414* * * 3415* * control_Nose_Te * 3416* * * 3417* **************************** 3418 real*8 function control_Nose_Te() 3419 implicit none 3420 3421* **** control_nose common block **** 3422 logical nose 3423 real*8 Pe,Te,Pr,Tr 3424 common / control_nblock / Pe,Te,Pr,Tr,nose 3425 3426 3427 control_Nose_Te = Te 3428 return 3429 end 3430 3431 3432* **************************** 3433* * * 3434* * control_Nose_Pr * 3435* * * 3436* **************************** 3437 real*8 function control_Nose_Pr() 3438 implicit none 3439 3440* **** control_nose common block **** 3441 logical nose 3442 real*8 Pe,Te,Pr,Tr 3443 common / control_nblock / Pe,Te,Pr,Tr,nose 3444 3445 3446 control_Nose_Pr = Pr 3447 return 3448 end 3449 3450 3451* **************************** 3452* * * 3453* * control_Nose_Tr * 3454* * * 3455* **************************** 3456 real*8 function control_Nose_Tr() 3457 implicit none 3458 3459* **** control_nose common block **** 3460 logical nose 3461 real*8 Pe,Te,Pr,Tr 3462 common / control_nblock / Pe,Te,Pr,Tr,nose 3463 3464 3465 control_Nose_Tr = Tr 3466 return 3467 end 3468 3469 3470* **************************** 3471* * * 3472* * control_Mulliken * 3473* * * 3474* **************************** 3475 logical function control_Mulliken() 3476 implicit none 3477 3478#include "bafdecls.fh" 3479#include "btdb.fh" 3480#include "errquit.fh" 3481#include "control.fh" 3482 3483* **** control_rtdb common block **** 3484 integer rtdb 3485 common / control_rtdb1 / rtdb 3486 3487* ***** local variables **** 3488 logical value 3489 3490 value = .false. 3491 if (code.eq.1) then 3492 if (.not.btdb_get(rtdb,'cpsd:mulliken',mt_log,1,value)) 3493 > value = .false. 3494 end if 3495 if ((code.eq.2).or.(code.eq.7)) then 3496 if (.not.btdb_get(rtdb,'cpmd:mulliken',mt_log,1,value)) 3497 > value = .false. 3498 end if 3499 3500 3501 if ((code.eq.3).or.(code.eq.11)) then 3502 if (.not.btdb_get(rtdb,'cgsd:mulliken',mt_log,1,value)) 3503 > value = .false. 3504 end if 3505 3506 if ((code.eq.5).or.(code.eq.11)) then 3507 if (.not.btdb_get(rtdb,'band:mulliken',mt_log,1,value)) 3508 > value = .false. 3509 end if 3510 3511 3512 control_Mulliken = value 3513 return 3514 end 3515 3516 3517* ***************************** 3518* * * 3519* * control_allow_translation * 3520* * * 3521* ***************************** 3522 logical function control_allow_translation() 3523 implicit none 3524 3525#include "bafdecls.fh" 3526#include "btdb.fh" 3527#include "inp.fh" 3528#include "util.fh" 3529 3530 logical value 3531 character*50 operation 3532 3533 logical control_qmmm 3534 external control_qmmm 3535 3536* **** control_rtdb common block **** 3537 integer rtdb 3538 common / control_rtdb1 / rtdb 3539 3540 if (.not.btdb_get(rtdb,'cgsd:allow_translation', 3541 > mt_log,1,value)) 3542 > value = .true. 3543 3544* *** read the current operation **** 3545 if (.not. btdb_cget(rtdb, 'task:operation', 1, operation)) 3546 $ operation = ' ' 3547 3548* *** allow translation of operation == freq||hessian *** 3549 if (inp_compare(.false.,'freq',operation)) value = .true. 3550 if (inp_compare(.false.,'hessian',operation)) value = .true. 3551 3552* *** allow translation for QMMM *** 3553 if (control_qmmm()) value = .true. 3554 3555 3556 control_allow_translation = value 3557 return 3558 end 3559 3560 3561 3562 3563* ***************************** 3564* * * 3565* * control_qmmm * 3566* * * 3567* ***************************** 3568 logical function control_qmmm() 3569 implicit none 3570 3571#include "bafdecls.fh" 3572#include "btdb.fh" 3573#include "util.fh" 3574 3575 logical task_qmmm 3576 3577* **** control_rtdb common block **** 3578 integer rtdb 3579 common / control_rtdb1 / rtdb 3580 3581 if( .not. btdb_get(rtdb,'task:QMMM',mt_log,1,task_qmmm)) 3582 > task_qmmm = .false. 3583 3584 control_qmmm = task_qmmm 3585 return 3586 end 3587 3588 3589* ***************************** 3590* * * 3591* * control_makehmass2 * 3592* * * 3593* ***************************** 3594 logical function control_makehmass2() 3595 implicit none 3596 3597#include "bafdecls.fh" 3598#include "btdb.fh" 3599#include "util.fh" 3600 3601 logical makehmass2 3602 3603* **** control_rtdb common block **** 3604 integer rtdb 3605 common / control_rtdb1 / rtdb 3606 3607 if( .not. btdb_get(rtdb,'nwpw:makehmass2',mt_log,1,makehmass2)) 3608 > makehmass2 = .true. 3609 3610 control_makehmass2 = makehmass2 3611 return 3612 end 3613 3614 3615 3616* ***************************** 3617* * * 3618* * control_MP2 * 3619* * * 3620* ***************************** 3621 logical function control_MP2() 3622 implicit none 3623 3624#include "bafdecls.fh" 3625#include "btdb.fh" 3626#include "util.fh" 3627 3628 logical mp2 3629 3630* **** control_rtdb common block **** 3631 integer rtdb 3632 common / control_rtdb1 / rtdb 3633 3634 if( .not. btdb_get(rtdb,'nwpw:MP2',mt_log,1,mp2)) 3635 > mp2 = .false. 3636 3637 control_MP2 = mp2 3638 return 3639 end 3640 3641 3642* ***************************** 3643* * * 3644* * control_2qintegrals * 3645* * * 3646* ***************************** 3647 logical function control_2qintegrals() 3648 implicit none 3649 3650#include "bafdecls.fh" 3651#include "btdb.fh" 3652#include "util.fh" 3653 3654 logical q2 3655 3656* **** control_rtdb common block **** 3657 integer rtdb 3658 common / control_rtdb1 / rtdb 3659 3660 if( .not. btdb_get(rtdb,'nwpw:2qintegrals',mt_log,1,q2)) 3661 > q2 = .false. 3662 3663 control_2qintegrals = q2 3664 return 3665 end 3666 3667 3668* **************************** 3669* * * 3670* * control_num_kvectors * 3671* * * 3672* **************************** 3673 integer function control_num_kvectors() 3674 implicit none 3675 3676#include "bafdecls.fh" 3677#include "btdb.fh" 3678#include "errquit.fh" 3679 3680* **** control_rtdb common block **** 3681 integer rtdb 3682 common / control_rtdb1 / rtdb 3683 3684* **** local variables **** 3685 logical value 3686 character*50 zone_name 3687 character*50 rtdb_name 3688 integer num_kvectors,l 3689 3690 if (.not.btdb_cget(rtdb,'band:zone_name',1,zone_name)) 3691 > zone_name = 'zone_default' 3692 3693 l = index(zone_name,' ') -1 3694 rtdb_name = zone_name(1:l)//':number_kvectors' 3695 if (.not.btdb_get(rtdb,rtdb_name,mt_int,1,num_kvectors)) 3696 > num_kvectors = 1 3697 3698 control_num_kvectors = num_kvectors 3699 return 3700 end 3701 3702* ********************************** 3703* * * 3704* * control_num_kvectors_print * 3705* * * 3706* ********************************** 3707 integer function control_num_kvectors_print() 3708 implicit none 3709 3710#include "bafdecls.fh" 3711#include "btdb.fh" 3712#include "errquit.fh" 3713 3714* **** control_rtdb common block **** 3715 integer rtdb 3716 common / control_rtdb1 / rtdb 3717 3718* **** local variables **** 3719 logical value 3720 character*50 zone_name 3721 character*50 rtdb_name 3722 integer num_kvectors_print,l 3723 3724 if (.not.btdb_cget(rtdb,'band:zone_name',1,zone_name)) 3725 > zone_name = 'zone_default' 3726 3727 l = index(zone_name,' ') -1 3728 rtdb_name = zone_name(1:l)//':number_kvectors_print' 3729 if (.not.btdb_get(rtdb,rtdb_name,mt_int,1,num_kvectors_print)) 3730 > num_kvectors_print = 0 3731 3732 control_num_kvectors_print = num_kvectors_print 3733 return 3734 end 3735 3736 3737* **************************** 3738* * * 3739* * control_ksvector * 3740* * * 3741* **************************** 3742 subroutine control_ksvector(i,ks) 3743 implicit none 3744 integer i 3745 real*8 ks(4) 3746 3747#include "bafdecls.fh" 3748#include "btdb.fh" 3749#include "errquit.fh" 3750 3751* **** control_rtdb common block **** 3752 integer rtdb 3753 common / control_rtdb1 / rtdb 3754 3755* **** local variables **** 3756 character*50 zone_name 3757 character*50 rtdb_name 3758 integer num_kvectors,l 3759 integer kvs(2) 3760 3761* **** external functions **** 3762 integer control_num_kvectors 3763 external control_num_kvectors 3764 3765 num_kvectors = control_num_kvectors() 3766 3767 if (.not.BA_push_get(mt_dbl,(4*num_kvectors),'kvs',kvs(2),kvs(1))) 3768 > call errquit('control_ksvector: out of stack', 0,MA_ERR) 3769 3770 if (.not.btdb_cget(rtdb,'band:zone_name',1,zone_name)) 3771 > zone_name = 'zone_default' 3772 3773 3774 l = index(zone_name,' ') -1 3775 rtdb_name = zone_name(1:l)//':kvectors' 3776 if (.not.btdb_get(rtdb,rtdb_name,mt_dbl, 3777 > (4*num_kvectors),dbl_mb(kvs(1)))) then 3778 call dcopy(4*num_kvectors,0.0d0,0,dbl_mb(kvs(1)),1) 3779 end if 3780 3781 ks(1) = dbl_mb(kvs(1)+4*(i-1)) 3782 ks(2) = dbl_mb(kvs(1)+4*(i-1)+1) 3783 ks(3) = dbl_mb(kvs(1)+4*(i-1)+2) 3784 ks(4) = dbl_mb(kvs(1)+4*(i-1)+3) 3785 3786 if (.not.BA_pop_stack(kvs(2))) 3787 > call errquit('control_ksvector: failed to free stack',0,MA_ERR) 3788 3789 return 3790 end 3791 3792* **************************** 3793* * * 3794* * control_kvector * 3795* * * 3796* **************************** 3797 subroutine control_kvector(i,kv) 3798 implicit none 3799 integer i 3800 real*8 kv(3) 3801 3802* **** local variables **** 3803 real*8 ks(4) 3804 3805* **** external functions **** 3806 real*8 lattice_unitg 3807 external lattice_unitg 3808 3809 call control_ksvector(i,ks) 3810 3811 kv(1) = ks(1)*lattice_unitg(1,1) 3812 > + ks(2)*lattice_unitg(1,2) 3813 > + ks(3)*lattice_unitg(1,3) 3814 kv(2) = ks(1)*lattice_unitg(2,1) 3815 > + ks(2)*lattice_unitg(2,2) 3816 > + ks(3)*lattice_unitg(2,3) 3817 kv(3) = ks(1)*lattice_unitg(3,1) 3818 > + ks(2)*lattice_unitg(3,2) 3819 > + ks(3)*lattice_unitg(3,3) 3820 3821 return 3822 end 3823 3824* ********************************** 3825* * * 3826* * control_monkhorst_pack_grid * 3827* * * 3828* ********************************** 3829 subroutine control_monkhorst_pack_grid(grid) 3830 implicit none 3831 integer grid(3) 3832 3833#include "bafdecls.fh" 3834#include "btdb.fh" 3835#include "errquit.fh" 3836 3837* **** control_rtdb common block **** 3838 integer rtdb 3839 common / control_rtdb1 / rtdb 3840 3841* **** local variables **** 3842 character*50 zone_name 3843 character*50 rtdb_name 3844 integer l 3845 3846 if (.not.btdb_cget(rtdb,'band:zone_name',1,zone_name)) 3847 > zone_name = 'zone_default' 3848 3849 l = index(zone_name,' ') -1 3850 rtdb_name = zone_name(1:l)//':monkhorst-pack' 3851 if (.not.btdb_get(rtdb,rtdb_name,mt_int,3,grid)) then 3852 grid(1) = 0 3853 grid(2) = 0 3854 grid(3) = 0 3855 end if 3856 3857 return 3858 end 3859 3860* ********************************** 3861* * * 3862* * control_fft_dos_grid * 3863* * * 3864* ********************************** 3865 subroutine control_fft_dos_grid(grid) 3866 implicit none 3867 integer grid(3) 3868 3869#include "bafdecls.fh" 3870#include "btdb.fh" 3871#include "errquit.fh" 3872 3873* **** control_rtdb common block **** 3874 integer rtdb 3875 common / control_rtdb1 / rtdb 3876 3877* **** local variables **** 3878 character*50 zone_name 3879 character*50 rtdb_name 3880 integer l 3881 3882 if (.not.btdb_cget(rtdb,'band_fft:zone_name', 3883 > 1,zone_name)) 3884 > zone_name = 'zone_fft_default' 3885 3886 l = index(zone_name,' ') -1 3887 rtdb_name = zone_name(1:l)//':dos-grid' 3888 if(.not.btdb_get(rtdb,rtdb_name,mt_int,3,grid)) then 3889 grid(1) = 0 3890 grid(2) = 0 3891 grid(3) = 0 3892 end if 3893 3894 return 3895 end 3896 3897 3898* ********************************** 3899* * * 3900* * control_ksvector_index * 3901* * * 3902* ********************************** 3903 integer function control_ksvector_index(ks) 3904 implicit none 3905 real*8 ks(*) 3906 3907#include "bafdecls.fh" 3908#include "btdb.fh" 3909#include "errquit.fh" 3910 3911* **** control_rtdb common block **** 3912 integer rtdb 3913 common / control_rtdb1 / rtdb 3914 3915* **** local variables **** 3916 character*50 zone_name 3917 character*50 rtdb_name 3918 integer num_kvectors,i,i1,i2 3919 integer kvs(2) 3920 real*8 d1,d2 3921 3922* **** external functions **** 3923 integer control_num_kvectors 3924 external control_num_kvectors 3925 3926 num_kvectors = control_num_kvectors() 3927 3928 if (.not.BA_push_get(mt_dbl,(4*num_kvectors),'kvs',kvs(2),kvs(1))) 3929 > call errquit('control_ksvector: out of stack', 0,MA_ERR) 3930 3931 if (.not.btdb_cget(rtdb,'band:zone_name',1,zone_name)) 3932 > zone_name = 'zone_default' 3933 3934 i = index(zone_name,' ') -1 3935 rtdb_name = zone_name(1:i)//':kvectors' 3936 if (.not.btdb_get(rtdb,rtdb_name,mt_dbl, 3937 > (4*num_kvectors),dbl_mb(kvs(1)))) then 3938 call dcopy(4*num_kvectors,0.0d0,0,dbl_mb(kvs(1)),1) 3939 end if 3940 3941 i1 = -1 3942 i2 = -1 3943 do i=1,num_kvectors 3944 d1 = dsqrt((dbl_mb(kvs(1)+4*(i-1)) - ks(1))**2 3945 > + (dbl_mb(kvs(1)+4*(i-1)+1) - ks(2))**2 3946 > + (dbl_mb(kvs(1)+4*(i-1)+2) - ks(3))**2) 3947 if (d1.lt.1.0e-6) i1 = i 3948 d2 = dsqrt((dbl_mb(kvs(1)+4*(i-1)) + ks(1))**2 3949 > + (dbl_mb(kvs(1)+4*(i-1)+1) + ks(2))**2 3950 > + (dbl_mb(kvs(1)+4*(i-1)+2) + ks(3))**2) 3951 if (d2.lt.1.0e-6) i2 = i 3952 end do 3953 3954 if (.not.BA_pop_stack(kvs(2))) 3955 > call errquit('control_ksvector: failed to free stack',0,MA_ERR) 3956 3957 if (i1.eq.-1) i1 = i2 3958 control_ksvector_index = i1 3959 return 3960 end 3961 3962 3963 3964 3965* ***************************** 3966* * * 3967* * control_TotalCharge * 3968* * * 3969* ***************************** 3970 real*8 function control_TotalCharge() 3971 implicit none 3972 3973#include "bafdecls.fh" 3974#include "btdb.fh" 3975 3976* **** control_rtdb common block **** 3977 integer rtdb 3978 common / control_rtdb1 / rtdb 3979 3980 double precision charge 3981 3982 charge = 0.0d0 3983 if (.not.btdb_get(rtdb,'charge',mt_dbl,1,charge)) then 3984 charge = 0.0d0 3985 end if 3986 3987 control_TotalCharge = charge 3988 return 3989 end 3990 3991 3992 3993* ************************************** 3994* * * 3995* * control_frac_occ_extra_orbitals * 3996* * * 3997* ************************************** 3998 integer function control_frac_occ_extra_orbitals() 3999 implicit none 4000 4001#include "bafdecls.fh" 4002#include "btdb.fh" 4003 4004* **** control_rtdb common block **** 4005 integer rtdb 4006 common / control_rtdb1 / rtdb 4007 4008 integer norbs 4009 4010 norbs = 0 4011 if (.not. 4012 > btdb_get(rtdb,'nwpw:frac_occ:extra_orbitals', 4013 > mt_int,1,norbs)) then 4014 norbs = 0 4015 end if 4016 4017 control_frac_occ_extra_orbitals = norbs 4018 return 4019 end 4020 4021 4022 4023 4024* ***************************** 4025* * * 4026* * control_rtdb * 4027* * * 4028* ***************************** 4029 integer function control_rtdb() 4030 implicit none 4031 4032* **** control_rtdb common block **** 4033 integer trtdb 4034 common / control_rtdb1 / trtdb 4035 4036 control_rtdb = trtdb 4037 return 4038 end 4039 4040* ***************************** 4041* * * 4042* * control_minimizer * 4043* * * 4044* ***************************** 4045 integer function control_minimizer() 4046 implicit none 4047 4048#include "control.fh" 4049 4050 control_minimizer = minimizer 4051 return 4052 end 4053 4054 4055 4056* ***************************** 4057* * * 4058* * control_lmbfgs_size * 4059* * * 4060* ***************************** 4061 integer function control_lmbfgs_size() 4062 implicit none 4063 4064#include "bafdecls.fh" 4065#include "btdb.fh" 4066 4067* **** control_rtdb common block **** 4068 integer rtdb 4069 common / control_rtdb1 / rtdb 4070 4071 integer lmbfgs_size 4072 4073 if (.not.btdb_get(rtdb,'nwpw:lmbfgs_size',mt_int,1,lmbfgs_size)) 4074 > then 4075 lmbfgs_size = 1 4076 end if 4077 4078 control_lmbfgs_size = lmbfgs_size 4079 return 4080 end 4081 4082* ***************************** 4083* * * 4084* * control_precondition * 4085* * * 4086* ***************************** 4087 logical function control_precondition() 4088 implicit none 4089 4090#include "control.fh" 4091 4092c#include "bafdecls.fh" 4093c#include "btdb.fh" 4094c 4095c* **** control_rtdb common block **** 4096c integer rtdb 4097c common / control_rtdb1 / rtdb 4098c 4099c logical precondition 4100c 4101c if (.not.btdb_get(rtdb,'nwpw:precondition', 4102c > mt_log,1,precondition)) 4103c > then 4104c precondition = .false. 4105c end if 4106 4107 control_precondition = precondition 4108 return 4109 end 4110 4111* ***************************** 4112* * * 4113* * control_lmbfgs_ondisk * 4114* * * 4115* ***************************** 4116 logical function control_lmbfgs_ondisk() 4117 implicit none 4118 4119#include "bafdecls.fh" 4120#include "btdb.fh" 4121 4122* **** control_rtdb common block **** 4123 integer rtdb 4124 common / control_rtdb1 / rtdb 4125 4126 logical ondisk 4127 4128 if (.not.btdb_get(rtdb,'nwpw:lmbfgs_ondisk', 4129 > mt_log,1,ondisk)) 4130 > then 4131 ondisk = .false. 4132 end if 4133 4134 control_lmbfgs_ondisk = ondisk 4135 return 4136 end 4137 4138* ***************************** 4139* * * 4140* * control_pspparameters * 4141* * * 4142* ***************************** 4143 subroutine control_pspparameters(atom,lmax,locp,rlocal) 4144 implicit none 4145 character*4 atom 4146 integer lmax,locp 4147 real*8 rlocal 4148 4149#include "bafdecls.fh" 4150#include "btdb.fh" 4151 4152* **** control_rtdb common block **** 4153 integer rtdb 4154 common / control_rtdb1 / rtdb 4155 4156 integer l,k 4157 character*5 element 4158 character*50 rtdb_name 4159 4160 element = ' ' 4161 element = atom 4162 l = index(element,' ') - 1 4163 rtdb_name = ' ' 4164 rtdb_name = element(1:l)//':lmax' 4165 k = index(rtdb_name,' ') - 1 4166 if (.not.btdb_get(rtdb,rtdb_name(1:k),mt_int,1,lmax)) 4167 > lmax = -1 4168 4169 rtdb_name = ' ' 4170 rtdb_name = element(1:l)//':locp' 4171 k = index(rtdb_name,' ') - 1 4172 if (.not.btdb_get(rtdb,rtdb_name(1:k),mt_int,1,locp)) 4173 > locp = -1 4174 4175 rtdb_name = ' ' 4176 rtdb_name = element(1:l)//':rlocal' 4177 k = index(rtdb_name,' ') - 1 4178 if (.not.btdb_get(rtdb,rtdb_name(1:k),mt_dbl,1,rlocal)) 4179 > rlocal = 1.0d0 4180 4181 return 4182 end 4183 4184* ********************************** 4185* * * 4186* * control_mullikenparameters * 4187* * * 4188* ********************************** 4189 subroutine control_mullikenparameters(atom,rcut,lmbda) 4190 implicit none 4191 character*4 atom 4192 real*8 rcut,lmbda 4193 4194#include "bafdecls.fh" 4195#include "btdb.fh" 4196 4197* **** control_rtdb common block **** 4198 integer rtdb 4199 common / control_rtdb1 / rtdb 4200 4201 logical value 4202 integer l,k 4203 character*5 element 4204 character*50 rtdb_name 4205 4206 element = ' ' 4207 element = atom 4208 l = index(element,' ') - 1 4209 rtdb_name = ' ' 4210 rtdb_name = element(1:l)//':mulliken:rcut' 4211 k = index(rtdb_name,' ') - 1 4212 if (.not.btdb_get(rtdb,rtdb_name(1:k),mt_dbl,1,rcut)) 4213 > rcut = 1.0d0 4214 4215 rtdb_name = ' ' 4216 rtdb_name = element(1:l)//':mulliken:lmbda' 4217 k = index(rtdb_name,' ') - 1 4218 if (.not.btdb_get(rtdb,rtdb_name(1:k),mt_dbl,1,lmbda)) 4219 > lmbda = 0.0d0 4220 4221 return 4222 end 4223 4224 4225 4226* *************************** 4227* * * 4228* * control_Ep * 4229* * * 4230* *************************** 4231 real*8 function control_Ep() 4232 implicit none 4233 4234#include "control.fh" 4235 4236 control_Ep = Ep 4237 return 4238 end 4239 4240 4241 4242* *************************** 4243* * * 4244* * control_Sp * 4245* * * 4246* *************************** 4247 real*8 function control_Sp() 4248 implicit none 4249 4250#include "control.fh" 4251 4252 control_Sp = Sp 4253 return 4254 end 4255 4256 4257* *************************** 4258* * * 4259* * control_SA * 4260* * * 4261* *************************** 4262 logical function control_SA() 4263 implicit none 4264 4265#include "control.fh" 4266 4267 control_SA=SA 4268 return 4269 end 4270 4271* *************************** 4272* * * 4273* * control_SA_decay * 4274* * * 4275* *************************** 4276 real*8 function control_SA_decay(choice) 4277 implicit none 4278 integer choice 4279 4280#include "control.fh" 4281 4282 control_SA_decay = sa_decay(choice) 4283 return 4284 end 4285 4286* ********************************** 4287* * * 4288* * control_dipole_motion * 4289* * * 4290* ********************************** 4291 logical function control_dipole_motion() 4292 implicit none 4293 4294#include "control.fh" 4295 4296 control_dipole_motion=dipole_motion 4297 return 4298 end 4299 4300* *************************** 4301* * * 4302* * control_Fei * 4303* * * 4304* *************************** 4305 logical function control_Fei() 4306 implicit none 4307 4308#include "control.fh" 4309 4310 control_Fei=fei 4311 return 4312 end 4313 4314 4315* *************************** 4316* * * 4317* * control_Fei_quench * 4318* * * 4319* *************************** 4320 logical function control_Fei_quench() 4321 implicit none 4322 4323#include "control.fh" 4324 4325 control_Fei_quench=fei_quench 4326 return 4327 end 4328 4329 4330* *************************** 4331* * * 4332* * control_gram_schmidt * 4333* * * 4334* *************************** 4335 logical function control_gram_schmidt() 4336 implicit none 4337 4338#include "control.fh" 4339 4340 control_gram_schmidt=gram_schmidt 4341 return 4342 end 4343 4344 4345* *************************** 4346* * * 4347* * control_balance * 4348* * * 4349* *************************** 4350 logical function control_balance() 4351 implicit none 4352 4353#include "control.fh" 4354 4355 control_balance=balance 4356 return 4357 end 4358 4359 4360 4361* ***************************** 4362* * * 4363* * control_excited_ne * 4364* * * 4365* ***************************** 4366 integer function control_excited_ne(ii) 4367 implicit none 4368 integer ii 4369 4370#include "bafdecls.fh" 4371#include "btdb.fh" 4372 4373* **** control_rtdb common block **** 4374 integer rtdb 4375 common / control_rtdb1 / rtdb 4376 4377 integer ne(2) 4378 4379 if (.not.btdb_get(rtdb,'nwpw:excited_ne',mt_int,2,ne)) then 4380 ne(1) = 0 4381 ne(2) = 0 4382 end if 4383 4384 control_excited_ne = ne(ii) 4385 return 4386 end 4387 4388 4389 4390* ************************************* 4391* * * 4392* * control_fractional_orbitals * 4393* * * 4394* ************************************* 4395 integer function control_fractional_orbitals(ii) 4396 implicit none 4397 integer ii 4398#include "control.fh" 4399 control_fractional_orbitals = frac_ne(ii) 4400 return 4401 end 4402* ************************************* 4403* * * 4404* * control_fractional_kT * 4405* * * 4406* ************************************* 4407 real*8 function control_fractional_kT() 4408 implicit none 4409#include "control.fh" 4410* *** local variables and parameters **** 4411 double precision kb 4412 parameter (kb=3.16679d-6) 4413 control_fractional_kT = kb*frac_temperature 4414 return 4415 end 4416* ************************************* 4417* * * 4418* * control_fractional_temperature * 4419* * * 4420* ************************************* 4421 real*8 function control_fractional_temperature() 4422 implicit none 4423#include "control.fh" 4424 control_fractional_temperature = frac_temperature 4425 return 4426 end 4427* ************************************* 4428* * * 4429* * control_fractional_smeartype * 4430* * * 4431* ************************************* 4432 integer function control_fractional_smeartype() 4433 implicit none 4434#include "control.fh" 4435 control_fractional_smeartype = frac_smeartype 4436 return 4437 end 4438* ************************************* 4439* * * 4440* * control_fractional * 4441* * * 4442* ************************************* 4443 logical function control_fractional() 4444 implicit none 4445#include "control.fh" 4446 control_fractional = fractional 4447 return 4448 end 4449 4450* ************************************* 4451* * * 4452* * control_ortho * 4453* * * 4454* ************************************* 4455 logical function control_ortho() 4456 implicit none 4457 4458#include "bafdecls.fh" 4459#include "btdb.fh" 4460 4461* **** control_rtdb common block **** 4462 integer rtdb 4463 common / control_rtdb1 / rtdb 4464 4465 logical ortho 4466 4467 if (.not. 4468 > btdb_get(rtdb,'nwpw:ortho_initialize',mt_log,1,ortho)) then 4469 ortho = .true. 4470 end if 4471 4472 control_ortho = ortho 4473 return 4474 end 4475 4476 4477 4478 4479* ***************************** 4480* * * 4481* * control_mapping * 4482* * * 4483* ***************************** 4484 integer function control_mapping() 4485 implicit none 4486 4487#include "control.fh" 4488 4489 control_mapping = mapping 4490 return 4491 end 4492 4493 4494* ***************************** 4495* * * 4496* * control_np_dimensions * 4497* * * 4498* ***************************** 4499 integer function control_np_dimensions(i) 4500 implicit none 4501 integer i 4502 4503#include "control.fh" 4504 4505 control_np_dimensions = np_dimensions(i) 4506 return 4507 end 4508 4509* ***************************** 4510* * * 4511* * control_np_orbital * 4512* * * 4513* ***************************** 4514 integer function control_np_orbital() 4515 implicit none 4516 4517#include "control.fh" 4518 4519 control_np_orbital = np_dimensions(2) 4520 return 4521 end 4522 4523 4524 4525* ***************************** 4526* * * 4527* * control_mapping1d * 4528* * * 4529* ***************************** 4530 integer function control_mapping1d() 4531 implicit none 4532 4533#include "control.fh" 4534 4535 control_mapping1d = mapping1d 4536 return 4537 end 4538 4539* ***************************** 4540* * * 4541* * control_parallel_io * 4542* * * 4543* ***************************** 4544 logical function control_parallel_io() 4545 implicit none 4546 4547#include "control.fh" 4548 4549 control_parallel_io = pio 4550 return 4551 end 4552 4553 4554 4555* ***************************** 4556* * * 4557* * control_oep * 4558* * * 4559* ***************************** 4560 logical function control_oep() 4561 implicit none 4562 4563#include "bafdecls.fh" 4564#include "btdb.fh" 4565 4566* **** control_rtdb common block **** 4567 integer rtdb 4568 common / control_rtdb1 / rtdb 4569 4570 logical oep 4571 4572 if (.not.btdb_get(rtdb,'nwpw:oep',mt_log,1,oep)) then 4573 oep = .false. 4574 end if 4575 4576 control_oep = oep 4577 return 4578 end 4579 4580 4581 4582* ***************************** 4583* * * 4584* * control_new_vpsi * 4585* * * 4586* ***************************** 4587 logical function control_new_vpsi() 4588 implicit none 4589 4590#include "bafdecls.fh" 4591#include "btdb.fh" 4592 4593* **** control_rtdb common block **** 4594 integer rtdb 4595 common / control_rtdb1 / rtdb 4596 4597 logical new_vpsi 4598 4599 if (.not.btdb_get(rtdb,'nwpw:new_vmovecs',mt_log,1,new_vpsi)) then 4600 new_vpsi = .false. 4601 end if 4602 4603 control_new_vpsi = new_vpsi 4604 return 4605 end 4606 4607 4608* ***************************** 4609* * * 4610* * control_COM_shift * 4611* * * 4612* ***************************** 4613 logical function control_COM_shift() 4614 implicit none 4615 4616#include "bafdecls.fh" 4617#include "btdb.fh" 4618 4619* **** control_rtdb common block **** 4620 integer rtdb 4621 common / control_rtdb1 / rtdb 4622 4623 logical com_shift 4624 4625 if (.not.btdb_get(rtdb,'nwpw:com_shift',mt_log,1,com_shift)) then 4626 com_shift = .true. 4627 end if 4628 4629 control_COM_shift = com_shift 4630 return 4631 end 4632 4633 4634 4635* ***************************** 4636* * * 4637* * control_DOS * 4638* * * 4639* ***************************** 4640 logical function control_DOS() 4641 implicit none 4642 4643#include "bafdecls.fh" 4644#include "btdb.fh" 4645 4646* **** control_rtdb common block **** 4647 integer rtdb 4648 common / control_rtdb1 / rtdb 4649 4650 logical dos 4651 real*8 alpha 4652 4653 dos = .false. 4654 if (btdb_get(rtdb,'dos:alpha',mt_dbl,1,alpha)) dos = .true. 4655 4656 control_DOS = dos 4657 return 4658 end 4659 4660 4661* ***************************** 4662* * * 4663* * control_psi_tmp * 4664* * * 4665* ***************************** 4666 logical function control_psi_tmp() 4667 implicit none 4668 4669#include "bafdecls.fh" 4670#include "btdb.fh" 4671 4672* **** control_rtdb common block **** 4673 integer rtdb 4674 common / control_rtdb1 / rtdb 4675 4676 logical psitmp 4677 4678 psitmp = .false. 4679 if (.not.btdb_get(rtdb,'nwpw:psi_tmp',mt_log,1,psitmp)) 4680 > psitmp = .false. 4681 4682 control_psi_tmp = psitmp 4683 return 4684 end 4685 4686* ***************************** 4687* * * 4688* * control_mulliken_kawai * 4689* * * 4690* ***************************** 4691 logical function control_mulliken_kawai() 4692 implicit none 4693 4694#include "bafdecls.fh" 4695#include "btdb.fh" 4696 4697* **** control_rtdb common block **** 4698 integer rtdb 4699 common / control_rtdb1 / rtdb 4700 4701 logical value 4702 4703 value = .false. 4704 if (.not.btdb_get(rtdb,'nwpw:mulliken_kawai',mt_log,1,value)) 4705 > value = .false. 4706 4707 control_mulliken_kawai = value 4708 return 4709 end 4710 4711 4712 4713* ***************************** 4714* * * 4715* * control_zero_forces * 4716* * * 4717* ***************************** 4718 logical function control_zero_forces() 4719 implicit none 4720 4721#include "bafdecls.fh" 4722#include "btdb.fh" 4723 4724* **** control_rtdb common block **** 4725 integer rtdb 4726 common / control_rtdb1 / rtdb 4727 4728 logical value 4729 4730 value = .false. 4731 if (.not.btdb_get(rtdb,'nwpw:zero_forces',mt_log,1,value)) 4732 > value = .false. 4733 4734 control_zero_forces = value 4735 return 4736 end 4737 4738 4739 4740* ******************************** 4741* * * 4742* * control_dos_grid_structure * 4743* * * 4744* ******************************** 4745 subroutine control_dos_grid_structure(grid) 4746 implicit none 4747 integer grid(3) 4748 4749#include "bafdecls.fh" 4750#include "btdb.fh" 4751 4752* **** control_rtdb common block **** 4753 integer rtdb 4754 common / control_rtdb1 / rtdb 4755 4756 if (.not.btdb_get(rtdb,'band:dos-grid',mt_int,3,grid)) then 4757 grid(1) = 1 4758 grid(2) = 1 4759 grid(3) = 1 4760 end if 4761 4762 return 4763 end 4764 4765 4766 4767 4768* *********************************** 4769* * * 4770* * control_reset_band_structure * 4771* * * 4772* *********************************** 4773 subroutine control_reset_band_structure() 4774 implicit none 4775 4776#include "bafdecls.fh" 4777#include "btdb.fh" 4778#include "control.fh" 4779 4780* **** control_rtdb common block **** 4781 integer rtdb 4782 common / control_rtdb1 / rtdb 4783 4784 if (.not. btdb_cget(rtdb, 'pspw:input bvectors', 4785 > 1,input_wavefunction_filename)) 4786 > input_wavefunction_filename = 'atomic' 4787 4788 if (.not. btdb_cget(rtdb, 'pspw:output bvectors', 4789 > 1,output_wavefunction_filename)) 4790 > output_wavefunction_filename = ' ' 4791 if (output_wavefunction_filename.eq.' ')then 4792 if (input_wavefunction_filename.eq.'atomic')then 4793 call util_file_prefix('bmovecs',output_wavefunction_filename) 4794 else 4795 output_wavefunction_filename = input_wavefunction_filename 4796 endif 4797 endif 4798 if (input_wavefunction_filename.eq.'atomic')then 4799 input_wavefunction_filename = output_wavefunction_filename 4800 end if 4801 4802 4803 return 4804 end 4805 4806* ************************************** 4807* * * 4808* * control_num_kvectors_structure * 4809* * * 4810* ************************************** 4811 integer function control_num_kvectors_structure() 4812 implicit none 4813 4814#include "bafdecls.fh" 4815#include "btdb.fh" 4816#include "errquit.fh" 4817 4818* **** control_rtdb common block **** 4819 integer rtdb 4820 common / control_rtdb1 / rtdb 4821 4822* **** local variables **** 4823 logical value 4824 character*50 zone_name 4825 character*50 rtdb_name 4826 integer num_kvectors,l 4827 4828 value = btdb_cget(rtdb,'band_structure:zone_name',1,zone_name) 4829 4830 l = index(zone_name,' ') -1 4831 rtdb_name = zone_name(1:l)//':number_kvectors' 4832 value = value.and. 4833 > btdb_get(rtdb,rtdb_name,mt_int,1,num_kvectors) 4834 4835 if (.not. value) 4836 > call errquit('control_num_kvectors_structure: failed', 4837 > 0, RTDB_ERR) 4838 4839 control_num_kvectors_structure = num_kvectors 4840 return 4841 end 4842 4843 4844* ************************************ 4845* * * 4846* * control_ksvector_structure * 4847* * * 4848* ************************************ 4849 subroutine control_ksvector_structure(i,ks) 4850 implicit none 4851 integer i 4852 real*8 ks(4) 4853 4854#include "bafdecls.fh" 4855#include "btdb.fh" 4856#include "errquit.fh" 4857 4858* **** control_rtdb common block **** 4859 integer rtdb 4860 common / control_rtdb1 / rtdb 4861 4862* **** local variables **** 4863 logical value 4864 character*50 zone_name 4865 character*50 rtdb_name 4866 integer num_kvectors,l 4867 integer kvs(2) 4868 4869* **** external functions **** 4870 integer control_num_kvectors_structure 4871 external control_num_kvectors_structure 4872 4873 num_kvectors = control_num_kvectors_structure() 4874 value = BA_push_get(mt_dbl,(4*num_kvectors), 4875 > 'kvs',kvs(2),kvs(1)) 4876 if (.not. value) 4877 > call errquit('control_ksvector: failed to get zone name', 0, 4878 & MA_ERR) 4879 4880 value = value.and. 4881 > btdb_cget(rtdb,'band_structure:zone_name',1,zone_name) 4882 if (.not. value) 4883 > call errquit('control_ksvector: failed to get zone name', 0, 4884 & RTDB_ERR) 4885 4886 l = index(zone_name,' ') -1 4887 rtdb_name = zone_name(1:l)//':kvectors' 4888 value = value.and. 4889 > btdb_get(rtdb,rtdb_name,mt_dbl, 4890 > (4*num_kvectors), 4891 > dbl_mb(kvs(1))) 4892 4893 if (.not. value) 4894 > call errquit('control_ksvector: failed to get kvs', 0, 4895 & RTDB_ERR) 4896 4897 ks(1) = dbl_mb(kvs(1)+4*(i-1)) 4898 ks(2) = dbl_mb(kvs(1)+4*(i-1)+1) 4899 ks(3) = dbl_mb(kvs(1)+4*(i-1)+2) 4900 ks(4) = dbl_mb(kvs(1)+4*(i-1)+3) 4901 4902 value = value.and.BA_pop_stack(kvs(2)) 4903 4904 if (.not. value) 4905 > call errquit('control_ksvector: failed to free stack', 0, 4906 & MA_ERR) 4907 return 4908 end 4909 4910 4911* ************************************ 4912* * * 4913* * control_kvector_structure * 4914* * * 4915* ************************************ 4916 subroutine control_kvector_structure(i,kv) 4917 implicit none 4918 integer i 4919 real*8 kv(3) 4920 4921* **** local variables **** 4922 real*8 ks(4) 4923 4924* **** external functions **** 4925 real*8 lattice_unitg 4926 external lattice_unitg 4927 4928 call control_ksvector_structure(i,ks) 4929 4930 kv(1) = ks(1)*lattice_unitg(1,1) 4931 > + ks(2)*lattice_unitg(1,2) 4932 > + ks(3)*lattice_unitg(1,3) 4933 kv(2) = ks(1)*lattice_unitg(2,1) 4934 > + ks(2)*lattice_unitg(2,2) 4935 > + ks(3)*lattice_unitg(2,3) 4936 kv(3) = ks(1)*lattice_unitg(3,1) 4937 > + ks(2)*lattice_unitg(3,2) 4938 > + ks(3)*lattice_unitg(3,3) 4939 4940 return 4941 end 4942 4943 4944 4945* ***************************** 4946* * * 4947* * control_print * 4948* * * 4949* ***************************** 4950 logical function control_print(level) 4951 implicit none 4952 integer level 4953 4954 4955* **** control_print common block **** 4956 integer print_level 4957 common / control_print1 / print_level 4958 4959 4960 logical value 4961 4962 if (level.le.print_level) then 4963 value = .true. 4964 else 4965 value = .false. 4966 end if 4967 4968 control_print = value 4969 return 4970 end 4971 4972 4973* ***************************** 4974* * * 4975* * control_reduce_print * 4976* * * 4977* ***************************** 4978 subroutine control_reduce_print() 4979 implicit none 4980 4981* **** control_print common block **** 4982 integer print_level 4983 common / control_print1 / print_level 4984 4985 print_level = print_level -10 4986 return 4987 end 4988 4989* ***************************** 4990* * * 4991* * control_up_print * 4992* * * 4993* ***************************** 4994 subroutine control_up_print() 4995 implicit none 4996 4997* **** control_print common block **** 4998 integer print_level 4999 common / control_print1 / print_level 5000 5001 print_level = print_level + 10 5002 return 5003 end 5004 5005 5006 5007 5008* ************************************ 5009* * * 5010* * control_optimize_cell_strategy * 5011* * * 5012* ************************************ 5013 5014 integer function control_optimize_cell_strategy() 5015 implicit none 5016 5017#include "bafdecls.fh" 5018#include "btdb.fh" 5019 5020* **** control_rtdb common block **** 5021 integer rtdb 5022 common / control_rtdb1 / rtdb 5023 5024 integer optimize_strategy 5025 5026 if (.not.btdb_get(rtdb,'cell_optimize:optimize_strategy', 5027 > mt_int,1,optimize_strategy)) 5028 > optimize_strategy = 0 5029 5030 control_optimize_cell_strategy = optimize_strategy 5031 return 5032 end 5033 5034 5035 5036 5037* ************************************* 5038* * * 5039* * control_optimize_lattice_vectors * 5040* * * 5041* ************************************* 5042 5043 integer function control_optimize_lattice_vectors(u,v) 5044 implicit none 5045 integer u,v 5046 5047#include "bafdecls.fh" 5048#include "btdb.fh" 5049 5050* **** control_rtdb common block **** 5051 integer rtdb 5052 common / control_rtdb1 / rtdb 5053 5054 integer optimize_lattice_vectors(3,3) 5055 5056 if (.not.btdb_get(rtdb,'cell_optimize:optimize_lattice_vectors', 5057 > mt_int,9,optimize_lattice_vectors)) then 5058 5059 optimize_lattice_vectors(1,1) =1 5060 optimize_lattice_vectors(2,1) =1 5061 optimize_lattice_vectors(3,1) =1 5062 optimize_lattice_vectors(1,2) =1 5063 optimize_lattice_vectors(2,2) =1 5064 optimize_lattice_vectors(3,2) =1 5065 optimize_lattice_vectors(1,3) =1 5066 optimize_lattice_vectors(2,3) =1 5067 optimize_lattice_vectors(3,3) =1 5068 end if 5069 5070 control_optimize_lattice_vectors = optimize_lattice_vectors(u,v) 5071 return 5072 end 5073 5074* ************************************* 5075* * * 5076* * control_optimize_lattice * 5077* * * 5078* ************************************* 5079 5080 integer function control_optimize_lattice(i) 5081 implicit none 5082 integer i 5083 5084#include "bafdecls.fh" 5085#include "btdb.fh" 5086 5087* **** control_rtdb common block **** 5088 integer rtdb 5089 common / control_rtdb1 / rtdb 5090 5091 integer optimize_lattice(6) 5092 5093 if (.not.btdb_get(rtdb,'cell_optimize:optimize_lattice', 5094 > mt_int,6,optimize_lattice)) then 5095 5096 optimize_lattice(1) =1 5097 optimize_lattice(2) =1 5098 optimize_lattice(3) =1 5099 optimize_lattice(4) =1 5100 optimize_lattice(5) =1 5101 optimize_lattice(6) =1 5102 end if 5103 5104 control_optimize_lattice = optimize_lattice(i) 5105 return 5106 end 5107 5108 5109 5110* ************************************ 5111* * * 5112* * control_lmax_multipole * 5113* * * 5114* ************************************ 5115 5116 integer function control_lmax_multipole() 5117 implicit none 5118 5119#include "bafdecls.fh" 5120#include "btdb.fh" 5121 5122* **** control_rtdb common block **** 5123 integer rtdb 5124 common / control_rtdb1 / rtdb 5125 5126 integer lmax 5127 5128 if (.not.btdb_get(rtdb,'nwpw:lmax_multipole',mt_int,1,lmax)) 5129 > lmax = 0 5130 5131 control_lmax_multipole = lmax 5132 return 5133 end 5134 5135 5136 5137* ************************************ 5138* * * 5139* * control_pfft3_qsize * 5140* * * 5141* ************************************ 5142 5143 integer function control_pfft3_qsize() 5144 implicit none 5145 5146#include "bafdecls.fh" 5147#include "btdb.fh" 5148 5149* **** control_rtdb common block **** 5150 integer rtdb 5151 common / control_rtdb1 / rtdb 5152 5153 integer qmax 5154 5155 if (.not.btdb_get(rtdb,'nwpw:pfft3_qsize',mt_int,1,qmax)) 5156 > qmax = 4 5157 5158 control_pfft3_qsize = qmax 5159 return 5160 end 5161 5162 5163* ************************************ 5164* * * 5165* * control_nprj_mult * 5166* * * 5167* ************************************ 5168 5169 integer function control_nprj_mult() 5170 implicit none 5171 5172#include "bafdecls.fh" 5173#include "btdb.fh" 5174 5175* **** control_rtdb common block **** 5176 integer rtdb 5177 common / control_rtdb1 / rtdb 5178 5179 integer qmax 5180 5181 if (.not.btdb_get(rtdb,'nwpw:nprj_mult',mt_int,1,qmax)) 5182 > qmax = 1 5183 5184 control_nprj_mult = qmax 5185 return 5186 end 5187 5188 5189 5190* ************************************ 5191* * * 5192* * control_symmetry * 5193* * * 5194* ************************************ 5195 5196 integer function control_symmetry() 5197 implicit none 5198 5199#include "control.fh" 5200 5201 control_symmetry = symm_number 5202 return 5203 end 5204 5205* ************************************ 5206* * * 5207* * control_spin_orbit * 5208* * * 5209* ************************************ 5210 5211 logical function control_spin_orbit() 5212 implicit none 5213 5214#include "control.fh" 5215 5216 control_spin_orbit = spin_orbit 5217 return 5218 end 5219 5220* ************************************ 5221* * * 5222* * control_fast_erf * 5223* * * 5224* ************************************ 5225 5226 logical function control_fast_erf() 5227 implicit none 5228 5229#include "control.fh" 5230 5231 control_fast_erf = fast_erf 5232 return 5233 end 5234 5235 5236 5237* ************************************ 5238* * * 5239* * control_fmm * 5240* * * 5241* ************************************ 5242 5243 logical function control_fmm() 5244 implicit none 5245 5246#include "control.fh" 5247 5248 control_fmm = fmm 5249 return 5250 end 5251 5252* ************************************ 5253* * * 5254* * control_periodic_dipole * 5255* * * 5256* ************************************ 5257 logical function control_periodic_dipole() 5258 implicit none 5259 5260#include "control.fh" 5261 5262 control_periodic_dipole = periodic_dipole 5263 return 5264 end 5265 5266* ************************************ 5267* * * 5268* * control_smooth_cutoff * 5269* * * 5270* ************************************ 5271 logical function control_smooth_cutoff() 5272 implicit none 5273 5274#include "control.fh" 5275 5276 control_smooth_cutoff = smooth_cutoff 5277 return 5278 end 5279 5280* ************************************ 5281* * * 5282* * control_smooth_cutoff_values * 5283* * * 5284* ************************************ 5285 real*8 function control_smooth_cutoff_values(i) 5286 implicit none 5287 integer i 5288 5289#include "control.fh" 5290 5291 control_smooth_cutoff_values = smooth_cutoff_values(i) 5292 return 5293 end 5294 5295 5296 5297 5298* ************************************ 5299* * * 5300* * control_fmm_lmax * 5301* * * 5302* ************************************ 5303 5304 integer function control_fmm_lmax() 5305 implicit none 5306 5307#include "control.fh" 5308 5309 control_fmm_lmax = fmm_lmax 5310 return 5311 end 5312 5313* ************************************ 5314* * * 5315* * control_fmm_lr * 5316* * * 5317* ************************************ 5318 5319 integer function control_fmm_lr() 5320 implicit none 5321 5322#include "control.fh" 5323 5324 control_fmm_lr = fmm_lr 5325 return 5326 end 5327 5328 5329 5330 5331 5332* ***************************** 5333* * * 5334* * control_pressure * 5335* * * 5336* ***************************** 5337 logical function control_pressure() 5338 implicit none 5339 5340#include "bafdecls.fh" 5341#include "btdb.fh" 5342 5343* **** control_rtdb common block **** 5344 integer rtdb 5345 common / control_rtdb1 / rtdb 5346 5347 logical value 5348 5349 value = .false. 5350 if (.not.btdb_get(rtdb,'cpmd:pressure',mt_log,1,value)) 5351 > value = .false. 5352 5353 control_pressure = value 5354 return 5355 end 5356 5357 5358* *********************************** 5359* * * 5360* * control_init_velocities * 5361* * * 5362* *********************************** 5363 logical function control_init_velocities() 5364 implicit none 5365 5366#include "bafdecls.fh" 5367#include "btdb.fh" 5368 5369* **** control_rtdb common block **** 5370 integer rtdb 5371 common / control_rtdb1 / rtdb 5372 5373 logical ok,tmp 5374 5375 if (.not.btdb_get(rtdb,'nwpw:init_velocities',mt_log,1,tmp)) 5376 > tmp = .false. 5377 5378 if (tmp) ok = rtdb_delete(rtdb,'nwpw:init_velocities') 5379 5380 control_init_velocities = tmp 5381 return 5382 end 5383 5384 5385* *********************************** 5386* * * 5387* * control_kbpp_ray * 5388* * * 5389* *********************************** 5390 logical function control_kbpp_ray() 5391 implicit none 5392 5393#include "bafdecls.fh" 5394#include "btdb.fh" 5395 5396* **** control_rtdb common block **** 5397 integer rtdb 5398 common / control_rtdb1 / rtdb 5399 5400 logical value 5401 5402 value = .true. 5403 if (.not.btdb_get(rtdb,'nwpw:kbpp_ray', 5404 > mt_log,1,value)) 5405 > value = .false. 5406 5407 control_kbpp_ray = value 5408 return 5409 end 5410 5411 5412 5413* *********************************** 5414* * * 5415* * control_kbpp_filter * 5416* * * 5417* *********************************** 5418 logical function control_kbpp_filter() 5419 implicit none 5420 5421#include "bafdecls.fh" 5422#include "btdb.fh" 5423 5424* **** control_rtdb common block **** 5425 integer rtdb 5426 common / control_rtdb1 / rtdb 5427 5428 logical value 5429 5430 value = .true. 5431 if (.not.btdb_get(rtdb,'nwpw:kbpp_filter', 5432 > mt_log,1,value)) 5433 > value = .false. 5434 5435 control_kbpp_filter = value 5436 return 5437 end 5438 5439 5440* *********************************** 5441* * * 5442* * control_brillioun_ondisk * 5443* * * 5444* *********************************** 5445 logical function control_brillioun_ondisk() 5446 implicit none 5447 5448#include "bafdecls.fh" 5449#include "btdb.fh" 5450 5451* **** control_rtdb common block **** 5452 integer rtdb 5453 common / control_rtdb1 / rtdb 5454 5455 logical value 5456 5457 value = .true. 5458 if (.not.btdb_get(rtdb,'nwpw:brillioun_ondisk', 5459 > mt_log,1,value)) 5460 > value = .false. 5461 5462 control_brillioun_ondisk = value 5463 return 5464 end 5465 5466 5467 5468* ************************************** 5469* * * 5470* * control_brillioun_use_symmetry * 5471* * * 5472* ************************************** 5473 logical function control_brillioun_use_symmetry() 5474 implicit none 5475 5476#include "bafdecls.fh" 5477#include "btdb.fh" 5478 5479* **** control_rtdb common block **** 5480 integer rtdb 5481 common / control_rtdb1 / rtdb 5482 5483 logical value 5484 5485 value = .true. 5486 if (.not.btdb_get(rtdb,'nwpw:brillioun_use_symmetry', 5487 > mt_log,1,value)) 5488 > value = .true. 5489 5490 control_brillioun_use_symmetry = value 5491 return 5492 end 5493 5494* ************************************** 5495* * * 5496* * control_rho_use_symmetry * 5497* * * 5498* ************************************** 5499 logical function control_rho_use_symmetry() 5500 implicit none 5501 5502#include "bafdecls.fh" 5503#include "btdb.fh" 5504 5505* **** control_rtdb common block **** 5506 integer rtdb 5507 common / control_rtdb1 / rtdb 5508 5509 logical value 5510 5511 value = .false. 5512 if (.not.btdb_get(rtdb,'nwpw:rho_use_symmetry', 5513 > mt_log,1,value)) 5514 > value = .false. 5515 5516 control_rho_use_symmetry = value 5517 return 5518 end 5519 5520 5521 5522 5523* *********************************** 5524* * * 5525* * control_mparallelized * 5526* * * 5527* *********************************** 5528 logical function control_mparallelized() 5529 implicit none 5530 5531#include "bafdecls.fh" 5532#include "btdb.fh" 5533 5534* **** control_rtdb common block **** 5535 integer rtdb 5536 common / control_rtdb1 / rtdb 5537 5538 logical value 5539 5540 value = .true. 5541 if (.not.btdb_get(rtdb,'nwpw:mparallelized',mt_log,1,value)) 5542 > value = .false. 5543 5544 control_mparallelized = value 5545 return 5546 end 5547 5548 5549* *********************************** 5550* * * 5551* * control_mreplicate_size * 5552* * * 5553* *********************************** 5554 integer function control_mreplicate_size() 5555 implicit none 5556 5557#include "bafdecls.fh" 5558#include "btdb.fh" 5559 5560* **** control_rtdb common block **** 5561 integer rtdb 5562 common / control_rtdb1 / rtdb 5563 5564 integer value 5565 5566 value = 16 5567 if (.not.btdb_get(rtdb,'nwpw:mreplicate_size',mt_int,1,value)) 5568 > value = 16 5569 5570 control_mreplicate_size = value 5571 return 5572 end 5573 5574 5575* *********************************** 5576* * * 5577* * control_bo_cpmd * 5578* * * 5579* *********************************** 5580 logical function control_bo_cpmd() 5581 implicit none 5582 5583#include "bafdecls.fh" 5584#include "btdb.fh" 5585 5586* **** control_rtdb common block **** 5587 integer rtdb 5588 common / control_rtdb1 / rtdb 5589 5590 logical value 5591 5592 value = .true. 5593 if (.not.btdb_get(rtdb,'nwpw:bo_cpmd',mt_log,1,value)) 5594 > value = .false. 5595 5596 control_bo_cpmd = value 5597 return 5598 end 5599 5600* *********************************** 5601* * * 5602* * control_hfxon_virtual * 5603* * * 5604* *********************************** 5605 logical function control_hfxon_virtual() 5606 implicit none 5607 5608#include "bafdecls.fh" 5609#include "btdb.fh" 5610 5611* **** control_rtdb common block **** 5612 integer rtdb 5613 common / control_rtdb1 / rtdb 5614 5615 logical value 5616 5617 value = .true. 5618 if (.not.btdb_get(rtdb,'nwpw:hfxon_virtual',mt_log,1,value)) 5619 > value = .true. 5620 5621 control_hfxon_virtual = value 5622 return 5623 end 5624 5625 5626* *********************************** 5627* * * 5628* * control_gradient_virtual * 5629* * * 5630* *********************************** 5631 logical function control_gradient_virtual() 5632 implicit none 5633 5634#include "bafdecls.fh" 5635#include "btdb.fh" 5636 5637* **** control_rtdb common block **** 5638 integer rtdb 5639 common / control_rtdb1 / rtdb 5640 5641 logical value 5642 5643 value = .false. 5644 if (.not.btdb_get(rtdb,'nwpw:gradient_virtual',mt_log,1,value)) 5645 > value = .false. 5646 5647 control_gradient_virtual = value 5648 return 5649 end 5650 5651* *********************************** 5652* * * 5653* * control_gradient_virtual_fac * 5654* * * 5655* *********************************** 5656 integer function control_gradient_virtual_fac() 5657 implicit none 5658 5659#include "bafdecls.fh" 5660#include "btdb.fh" 5661 5662* **** control_rtdb common block **** 5663 integer rtdb 5664 common / control_rtdb1 / rtdb 5665 5666 integer fac 5667 5668 fac = 19 5669 if (.not.btdb_get(rtdb,'nwpw:gradient_virtual_fac',mt_int,1,fac)) 5670 > fac = 19 5671 5672 control_gradient_virtual_fac = fac 5673 return 5674 end 5675 5676 5677* *********************************** 5678* * * 5679* * control_bound_virtual * 5680* * * 5681* *********************************** 5682 logical function control_bound_virtual() 5683 implicit none 5684 5685#include "bafdecls.fh" 5686#include "btdb.fh" 5687 5688* **** control_rtdb common block **** 5689 integer rtdb 5690 common / control_rtdb1 / rtdb 5691 5692 logical value 5693 5694 value = .false. 5695 if (.not.btdb_get(rtdb,'nwpw:bound_virtual',mt_log,1,value)) 5696 > value = .false. 5697 5698 control_bound_virtual = value 5699 return 5700 end 5701 5702* *********************************** 5703* * * 5704* * control_2x2_virtual * 5705* * * 5706* *********************************** 5707 logical function control_2x2_virtual() 5708 implicit none 5709 5710#include "bafdecls.fh" 5711#include "btdb.fh" 5712 5713* **** control_rtdb common block **** 5714 integer rtdb 5715 common / control_rtdb1 / rtdb 5716 5717 logical value 5718 5719 value = .false. 5720 if (.not.btdb_get(rtdb,'nwpw:2x2_virtual',mt_log,1,value)) 5721 > value = .false. 5722 5723 control_2x2_virtual = value 5724 return 5725 end 5726 5727* *********************************** 5728* * * 5729* * control_2x2ne_virtual * 5730* * * 5731* *********************************** 5732 logical function control_2x2ne_virtual() 5733 implicit none 5734 5735#include "bafdecls.fh" 5736#include "btdb.fh" 5737 5738* **** control_rtdb common block **** 5739 integer rtdb 5740 common / control_rtdb1 / rtdb 5741 5742 logical value 5743 5744 value = .false. 5745 if (.not.btdb_get(rtdb,'nwpw:2x2ne_virtual',mt_log,1,value)) 5746 > value = .false. 5747 5748 control_2x2ne_virtual = value 5749 return 5750 end 5751 5752 5753 5754* *********************************** 5755* * * 5756* * control_3x3_virtual * 5757* * * 5758* *********************************** 5759 logical function control_3x3_virtual() 5760 implicit none 5761 5762#include "bafdecls.fh" 5763#include "btdb.fh" 5764 5765* **** control_rtdb common block **** 5766 integer rtdb 5767 common / control_rtdb1 / rtdb 5768 5769 logical value 5770 5771 value = .false. 5772 if (.not.btdb_get(rtdb,'nwpw:3x3_virtual',mt_log,1,value)) 5773 > value = .false. 5774 5775 control_3x3_virtual = value 5776 return 5777 end 5778 5779* *********************************** 5780* * * 5781* * control_4x4_virtual * 5782* * * 5783* *********************************** 5784 logical function control_4x4_virtual() 5785 implicit none 5786 5787#include "bafdecls.fh" 5788#include "btdb.fh" 5789 5790* **** control_rtdb common block **** 5791 integer rtdb 5792 common / control_rtdb1 / rtdb 5793 5794 logical value 5795 5796 value = .false. 5797 if (.not.btdb_get(rtdb,'nwpw:4x4_virtual',mt_log,1,value)) 5798 > value = .false. 5799 5800 control_4x4_virtual = value 5801 return 5802 end 5803 5804 5805* *********************************** 5806* * * 5807* * control_CI_virtual * 5808* * * 5809* *********************************** 5810 integer function control_CI_virtual() 5811 implicit none 5812 5813#include "bafdecls.fh" 5814#include "btdb.fh" 5815 5816* **** control_rtdb common block **** 5817 integer rtdb 5818 common / control_rtdb1 / rtdb 5819 5820 integer ivalue 5821 5822 logical control_2x2_virtual,control_3x3_virtual 5823 external control_2x2_virtual,control_3x3_virtual 5824 logical control_4x4_virtual 5825 external control_4x4_virtual 5826 logical control_2x2ne_virtual 5827 external control_2x2ne_virtual 5828 5829 ivalue = 0 5830 if (.not.btdb_get(rtdb,'nwpw:CI_virtual',mt_int,1,ivalue)) then 5831 if (control_2x2_virtual()) then 5832 ivalue = 2 5833 else if (control_3x3_virtual()) then 5834 ivalue = 3 5835 else if (control_4x4_virtual()) then 5836 ivalue = 4 5837 else if (control_2x2ne_virtual()) then 5838 ivalue = 5 5839 end if 5840 end if 5841 5842 control_CI_virtual = ivalue 5843 return 5844 end 5845 5846 5847 5848 5849 5850* *********************************** 5851* * * 5852* * control_confine_virtual * 5853* * * 5854* *********************************** 5855 logical function control_confine_virtual() 5856 implicit none 5857 5858#include "bafdecls.fh" 5859#include "btdb.fh" 5860 5861* **** control_rtdb common block **** 5862 integer rtdb 5863 common / control_rtdb1 / rtdb 5864 5865 logical value 5866 5867 value = .false. 5868 if (.not.btdb_get(rtdb,'nwpw:confine_virtual',mt_log,1,value)) 5869 > value = .false. 5870 5871 control_confine_virtual = value 5872 return 5873 end 5874 5875 5876 5877 5878 5879 5880 5881 5882* ********************************************** 5883* * * 5884* * control_init_velocities_temperature * 5885* * * 5886* ********************************************** 5887 real*8 function control_init_velocities_temperature() 5888 implicit none 5889 5890#include "bafdecls.fh" 5891#include "btdb.fh" 5892 5893* **** control_rtdb common block **** 5894 integer rtdb 5895 common / control_rtdb1 / rtdb 5896 5897 real*8 value 5898 5899 if (.not.btdb_get(rtdb,'nwpw:init_velocities_temperature', 5900 > mt_dbl,1,value)) 5901 > value = 300.0d0 5902 5903 control_init_velocities_temperature = value 5904 return 5905 end 5906 5907* ********************************************** 5908* * * 5909* * control_init_velocities_seed * 5910* * * 5911* ********************************************** 5912 integer function control_init_velocities_seed() 5913 implicit none 5914 5915#include "bafdecls.fh" 5916#include "btdb.fh" 5917 5918* **** control_rtdb common block **** 5919 integer rtdb 5920 common / control_rtdb1 / rtdb 5921 5922 integer seed 5923 5924 if (.not.btdb_get(rtdb,'nwpw:init_velocities_seed', 5925 > mt_int,1,seed)) 5926 > seed = 494 5927 5928 control_init_velocities_seed = seed 5929 return 5930 end 5931 5932 5933 5934 5935* *********************************** 5936* * * 5937* * control_wannier_timestep * 5938* * * 5939* *********************************** 5940 real*8 function control_wannier_timestep() 5941 implicit none 5942 5943#include "bafdecls.fh" 5944#include "btdb.fh" 5945 5946* **** control_rtdb common block **** 5947 integer rtdb 5948 common / control_rtdb1 / rtdb 5949 5950 real*8 value 5951 5952 if (.not.btdb_get(rtdb,'wannier:time_step',mt_dbl,1,value)) 5953 > value = 2.7e-2 5954 5955 control_wannier_timestep = value 5956 return 5957 end 5958 5959 5960 5961* *********************************** 5962* * * 5963* * control_wannier_maxiter * 5964* * * 5965* *********************************** 5966 integer function control_wannier_maxiter() 5967 implicit none 5968 5969#include "bafdecls.fh" 5970#include "btdb.fh" 5971 5972* **** control_rtdb common block **** 5973 integer rtdb 5974 common / control_rtdb1 / rtdb 5975 5976 integer value 5977 5978 if (.not.btdb_get(rtdb,'wannier:maxiter',mt_int,1,value)) 5979 > value = 500 5980 5981 control_wannier_maxiter = value 5982 return 5983 end 5984 5985 5986* ********************************************** 5987* * * 5988* * control_attenuation * 5989* * * 5990* ********************************************** 5991 real*8 function control_attenuation() 5992 implicit none 5993#include "control.fh" 5994 control_attenuation = attenuation 5995 return 5996 end 5997 5998 5999 6000ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 6001 subroutine set_two_component_pseudopotential() 6002ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 6003c this called to signal that a two_component ppot is in use 6004ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 6005#include "control.fh" 6006 two_comp_ppot=.true. 6007 return 6008 end 6009ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 6010 logical function two_component_pseudopotential() 6011#include "control.fh" 6012 two_component_pseudopotential=two_comp_ppot 6013 return 6014 end 6015ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 6016 6017 6018 6019 6020* *********************************** 6021* * * 6022* * control_ecut_wcut_default * 6023* * * 6024* *********************************** 6025 subroutine control_ecut_wcut_default(rtdb,ecut,wcut) 6026 implicit none 6027 integer rtdb 6028 real*8 ecut,wcut 6029 6030#include "bafdecls.fh" 6031#include "btdb.fh" 6032#include "beom.fh" 6033#include "errquit.fh" 6034 6035* **** local variables **** 6036 integer taskid,MASTER 6037 parameter (MASTER=0) 6038 6039 logical value 6040 integer geom,ii,nion,l,h 6041 real*8 wcut_max,wcut_atom,q,rion(3) 6042 character*4 aname 6043 character*16 t,rtdbname 6044 character*255 sdir_name 6045 6046 6047* **** external functions **** 6048 real*8 nwpw_libgcutoff 6049 external nwpw_libgcutoff 6050 character*4 ion_aname_geom 6051 external ion_aname_geom 6052 logical parsepointcharge 6053 external parsepointcharge 6054 6055 !*** parse psp information on rtdb *** 6056 if (.not.btdb_get(rtdb,'nwpw:psp:cutoff',mt_dbl,1,wcut_max)) then 6057 call Parallel_taskid(taskid) 6058 wcut_max = -1.0d0 6059 value = beom_create(geom,'geometry') 6060 value = value.and.beom_rtdb_load(rtdb,geom,'geometry') 6061 value = value.and.geom_ncent(geom,nion) 6062 if (.not. value) 6063 > call errquit('control_ecut_wcut_default:cannot load geometry', 6064 > 0,GEOM_ERR) 6065 6066 do ii=1,nion 6067 if (.not.geom_cent_get(geom,ii,t,rion,q)) 6068 > call errquit('control_ecut_wcut_default:error reading ions', 6069 > 0,GEOM_ERR) 6070 if (.not.parsepointcharge(t)) then 6071 aname = ion_aname_geom(geom,ii) 6072 l = index(aname,' ') - 1 6073 if (l.le.0) l = 4 6074 rtdbname = aname(1:l)//':cutoff' 6075 l = l+7 6076 if (.not.btdb_get(rtdb,rtdbname(1:l),mt_dbl,1,wcut_atom)) 6077 > then 6078 6079 !*** define wcut_atom from psplibrary *** 6080 value = btdb_parallel(.false.) 6081 if (taskid.eq.MASTER) then 6082 call util_directory_name(sdir_name,.true.,0) 6083 h = index(sdir_name,' ') - 1 6084 open(unit=99,file=sdir_name(1:h)//'/junk.inp', 6085 > status='unknown') 6086 close(unit=99,status='delete') 6087 6088 call nwpw_libgeninp(1,aname, 6089 > sdir_name(1:h)//'/junk.inp') 6090 wcut_atom=nwpw_libgcutoff(aname) 6091 6092 if(.not.btdb_put(rtdb,rtdbname(1:l), 6093 > mt_dbl,1,wcut_atom)) 6094 > call errquit( 6095 > 'control_ecut_wcut_default:cannot write wcut_atom',0,0) 6096 6097 end if 6098 value = btdb_parallel(.true.) 6099 call Parallel_Brdcst_value(MASTER,wcut_atom) 6100 6101 end if 6102 6103 !*** reset wcut_max *** 6104 if (wcut_atom.gt.wcut_max) then 6105 wcut_max = wcut_atom 6106 value = btdb_parallel(.false.) 6107 if (taskid.eq.MASTER) then 6108 if (.not.btdb_put(rtdb,'nwpw:psp:cutoff', 6109 > mt_dbl,1,wcut_max)) 6110 > call errquit( 6111 > 'control_ecut_wcut_default:cannot write wcut_max',0,0) 6112 end if 6113 value = btdb_parallel(.true.) 6114 end if 6115 end if 6116 end do 6117 if (.not. beom_destroy(geom)) 6118 > call errquit('control_ecut_wcut_default:cannot destroy geom', 6119 > 0,GEOM_ERR) 6120 end if 6121 6122 wcut = wcut_max 6123 ecut = 2.0d0*wcut 6124 6125 return 6126 end 6127 6128 6129 6130* *********************************** 6131* * * 6132* * control_pspspin * 6133* * * 6134* *********************************** 6135 logical function control_pspspin() 6136 implicit none 6137 6138#include "bafdecls.fh" 6139#include "btdb.fh" 6140 6141* **** control_rtdb common block **** 6142 integer rtdb 6143 common / control_rtdb1 / rtdb 6144 6145 logical value 6146 6147 value = .true. 6148 if (.not.btdb_get(rtdb,'nwpw:pspspin',mt_log,1,value)) 6149 > value = .false. 6150 6151 control_pspspin = value 6152 return 6153 end 6154 6155 6156* *********************************** 6157* * * 6158* * control_set_pspspin * 6159* * * 6160* *********************************** 6161 subroutine control_set_pspspin(nion, 6162 > upscale,downscale, 6163 > upl,downl, 6164 > upm,downm, 6165 > upions,downions) 6166 implicit none 6167 integer nion 6168 real*8 upscale(*),downscale(*) 6169 integer upl(*),downl(*),upm(*),downm(*) 6170 logical upions(*),downions(*) 6171 6172#include "bafdecls.fh" 6173#include "btdb.fh" 6174#include "errquit.fh" 6175 6176* **** control_rtdb common block **** 6177 integer rtdb 6178 common / control_rtdb1 / rtdb 6179 6180 integer taskid,MASTER 6181 parameter (MASTER=0) 6182 logical iamup 6183 integer i,ma_type,nactive_atoms,h_actlist,l_actlist 6184 integer pcount,ip,tl,tm 6185 real*8 tscale 6186 character*50 rtdb_name 6187 6188* *** external functions *** 6189 character*7 c_index_name 6190 external c_index_name 6191 6192 call Parallel_taskid(taskid) 6193 if (.not.btdb_get(rtdb,'nwpw:pspspin_count',mt_int,1,pcount)) 6194 > pcount = 0 6195 6196 if (taskid.eq.MASTER) write(*,2300) 6197 do ip=1,pcount 6198 6199 rtdb_name = 'nwpw:pspspin_iamup:'//c_index_name(ip) 6200 if (.not.btdb_get(rtdb,rtdb_name,mt_log,1,iamup)) iamup=.true. 6201 6202 if (iamup) then 6203 rtdb_name = 'nwpw:pspspin_upscale:'//c_index_name(ip) 6204 if (.not.btdb_get(rtdb,rtdb_name,mt_dbl,1,tscale)) 6205 > tscale = 1.0d0 6206 rtdb_name = 'nwpw:pspspin_upl:'//c_index_name(ip) 6207 if (.not.btdb_get(rtdb,rtdb_name,mt_int,1,tl)) 6208 > tl = -1 6209 rtdb_name = 'nwpw:pspspin_upm:'//c_index_name(ip) 6210 if (.not.btdb_get(rtdb,rtdb_name,mt_int,1,tm)) 6211 > tm = 99999 6212 if (taskid.eq.MASTER) then 6213 if (tm.lt.999) then 6214 write(*,2305) tl,tm,tscale 6215 else 6216 write(*,2301) tl,tscale 6217 end if 6218 end if 6219 6220 rtdb_name = 'nwpw:pspspin_upions:'//c_index_name(ip) 6221 if (rtdb_ma_get(rtdb,rtdb_name, ma_type, 6222 > nactive_atoms, h_actlist)) then 6223 if (.not.BA_get_index(h_actlist,l_actlist)) 6224 > call errquit( 6225 > 'control_set_pspspin: ma_get_index failed',0, 6226 > MA_ERR) 6227 6228 if (taskid.eq.MASTER) 6229 > write(*,2302) (int_mb(l_actlist+i),i=0,nactive_atoms-1) 6230 6231 do i=1,nactive_atoms 6232 upions(int_mb(l_actlist+i-1)) = .true. 6233 upl(int_mb(l_actlist+i-1)) = tl 6234 upm(int_mb(l_actlist+i-1)) = tm 6235 upscale(int_mb(l_actlist+i-1)) = tscale 6236 end do 6237 if (.not. BA_free_heap(h_actlist)) 6238 > call errquit( 6239 > 'control_set_pspspin:error freeing heap memory',0,MA_ERR) 6240 end if 6241 6242 6243 else 6244 rtdb_name = 'nwpw:pspspin_downscale:'//c_index_name(ip) 6245 if (.not.btdb_get(rtdb,rtdb_name,mt_dbl,1,tscale)) 6246 > tscale = 1.0d0 6247 rtdb_name = 'nwpw:pspspin_downl:'//c_index_name(ip) 6248 if (.not.btdb_get(rtdb,rtdb_name,mt_int,1,tl)) 6249 > tl = -1 6250 rtdb_name = 'nwpw:pspspin_downm:'//c_index_name(ip) 6251 if (.not.btdb_get(rtdb,rtdb_name,mt_int,1,tm)) 6252 > tm = 99999 6253 if (taskid.eq.MASTER) then 6254 if (tm.lt.999) then 6255 write(*,2306) tl,tm,tscale 6256 else 6257 write(*,2303) tl,tscale 6258 end if 6259 end if 6260 6261 rtdb_name = 'nwpw:pspspin_downions:'//c_index_name(ip) 6262 if (rtdb_ma_get(rtdb,rtdb_name,ma_type, 6263 > nactive_atoms,h_actlist)) then 6264 if (.not.BA_get_index(h_actlist,l_actlist)) 6265 > call errquit( 6266 > 'control_set_pspspin: ma_get_index failed',1, 6267 > MA_ERR) 6268 6269 if (taskid.eq.MASTER) 6270 > write(*,2304) (int_mb(l_actlist+i),i=0,nactive_atoms-1) 6271 6272 do i=1,nactive_atoms 6273 downions(int_mb(l_actlist+i-1)) = .true. 6274 downl(int_mb(l_actlist+i-1)) = tl 6275 downm(int_mb(l_actlist+i-1)) = tm 6276 downscale(int_mb(l_actlist+i-1)) = tscale 6277 end do 6278 if (.not.BA_free_heap(h_actlist)) 6279 > call errquit( 6280 > 'control_set_pspspin:error freeing heap memory',1,MA_ERR) 6281 end if 6282 end if 6283 end do 6284 6285 2300 format(/1x,"Antiferromagnetic Pentalty Function Input:") 6286 2301 format(2x," - pspspin: up l =",I2,10x," scale =",F8.3) 6287 2302 format(2x," - pspspin: up ion indexes =",10I5) 6288 2303 format(2x," - pspspin: down l =",I2,10x," scale =",F8.3) 6289 2304 format(2x," - pspspin: down ion indexes =",10I5) 6290 2305 format(2x," - pspspin: up l =",I2," not_m=",I3," scale =",F8.3) 6291 2306 format(2x," - pspspin: down l =",I2," not_m=",I3," scale =",F8.3) 6292 6293 return 6294 end 6295 6296* *********************************** 6297* * * 6298* * control_psputerm * 6299* * * 6300* *********************************** 6301 logical function control_psputerm() 6302 implicit none 6303 6304#include "bafdecls.fh" 6305#include "btdb.fh" 6306 6307* **** control_rtdb common block **** 6308 integer rtdb 6309 common / control_rtdb1 / rtdb 6310 6311 logical value 6312 6313 value = .true. 6314 if (.not.btdb_get(rtdb,'nwpw:uterm',mt_log,1,value)) 6315 > value = .false. 6316 6317 control_psputerm = value 6318 return 6319 end 6320 6321* *********************************** 6322* * * 6323* * control_pspnuterms * 6324* * * 6325* *********************************** 6326 integer function control_pspnuterms() 6327 implicit none 6328 6329#include "bafdecls.fh" 6330#include "btdb.fh" 6331 6332* **** control_rtdb common block **** 6333 integer rtdb 6334 common / control_rtdb1 / rtdb 6335 6336 integer nuterms 6337 6338 nuterms = 0 6339 if (.not.btdb_get(rtdb,'nwpw:nuterms',mt_int,1,nuterms)) 6340 > nuterms = 0 6341 6342 control_pspnuterms = nuterms 6343 return 6344 end 6345 6346 6347* *********************************** 6348* * * 6349* * control_set_psputerm * 6350* * * 6351* *********************************** 6352 subroutine control_set_psputerm(nion,nuterms,l,uscale,jscale,ions) 6353 implicit none 6354 integer nion,nuterms,l(nuterms) 6355 real*8 uscale(nuterms),jscale(nuterms) 6356 logical ions(nion,nuterms) 6357 6358#include "bafdecls.fh" 6359#include "btdb.fh" 6360#include "errquit.fh" 6361 6362* **** control_rtdb common block **** 6363 integer rtdb 6364 common / control_rtdb1 / rtdb 6365 6366 integer taskid,MASTER 6367 parameter (MASTER=0) 6368 integer nu,i,ma_type,nactive_atoms,h_actlist,l_actlist 6369 character*50 rtdb_name 6370 6371* **** external functions **** 6372 character*7 c_index_name 6373 external c_index_name 6374 6375 call Parallel_taskid(taskid) 6376 6377 if (taskid.eq.MASTER) write(*,3300) 6378 6379 do nu=1,nuterms 6380 6381 rtdb_name = 'nwpw:uterm_scale:'//c_index_name(nu) 6382 if (.not.btdb_get(rtdb,rtdb_name,mt_dbl,1,uscale(nu))) 6383 > uscale(nu) = 0.0d0 6384 6385 rtdb_name = 'nwpw:jterm_scale:'//c_index_name(nu) 6386 if (.not.btdb_get(rtdb,rtdb_name,mt_dbl,1,jscale(nu))) 6387 > jscale(nu) = 0.0d0 6388 6389 rtdb_name = 'nwpw:uterm_l:'//c_index_name(nu) 6390 if (.not.btdb_get(rtdb,rtdb_name,mt_int,1,l(nu))) 6391 > l(nu) = -1 6392 6393 do i=1,nion 6394 ions(i,nu) = .false. 6395 end do 6396 6397 if (taskid.eq.MASTER) write(*,3301) l(nu),uscale(nu),jscale(nu) 6398 6399 rtdb_name = 'nwpw:uterm_ions:'//c_index_name(nu) 6400 if (rtdb_ma_get(rtdb, rtdb_name,ma_type, 6401 > nactive_atoms, h_actlist)) then 6402 if (.not.BA_get_index(h_actlist,l_actlist)) 6403 > call errquit( 6404 > 'control_set_psputerm: ma_get_index failed',0, 6405 > MA_ERR) 6406 6407 if (taskid.eq.MASTER) then 6408 write(*,3302) (int_mb(l_actlist+i),i=0,nactive_atoms-1) 6409 write(*,*) 6410 end if 6411 6412 do i=1,nactive_atoms 6413 ions(int_mb(l_actlist+i-1),nu) = .true. 6414 end do 6415 if (.not. BA_free_heap(h_actlist)) 6416 > call errquit( 6417 > 'control_set_psputerm:error freeing heap memory',0,MA_ERR) 6418 end if 6419 6420 end do 6421 6422 3300 format(/1x,"Hubbard Uterm Function Input:") 6423 3301 format(2x," - uterm: l =",I2," U=",F8.3," J=",F8.3) 6424 3302 format(2x," - uterm: ion indexes =",10I5) 6425 return 6426 end 6427 6428 6429 6430* *********************************** 6431* * * 6432* * control_use_grid_cmp * 6433* * * 6434* *********************************** 6435 logical function control_use_grid_cmp() 6436 implicit none 6437 6438#include "bafdecls.fh" 6439#include "btdb.fh" 6440 6441* **** control_rtdb common block **** 6442 integer rtdb 6443 common / control_rtdb1 / rtdb 6444 6445 logical value 6446 6447 value = .true. 6448 if (.not.btdb_get(rtdb,'nwpw:use_grid_cmp',mt_log,1,value)) 6449 > value = .false. 6450 6451 control_use_grid_cmp = value 6452 return 6453 end 6454 6455 6456* *********************************** 6457* * * 6458* * control_wgc_alphabetalambda * 6459* * * 6460* *********************************** 6461 real*8 function control_wgc_alphabetalambda(i) 6462 implicit none 6463 integer i 6464 6465#include "bafdecls.fh" 6466#include "btdb.fh" 6467 6468* **** control_rtdb common block **** 6469 integer rtdb 6470 common / control_rtdb1 / rtdb 6471 6472 real*8 value 6473 6474 if (i.eq.1) then 6475 if (.not.btdb_get(rtdb,'nwpw:wgc_alpha',mt_dbl,1,value)) 6476 > value = 5.0d0/6.0d0 6477 else if (i.eq.2) then 6478 if (.not.btdb_get(rtdb,'nwpw:wgc_beta',mt_dbl,1,value)) 6479 > value = 5.0d0/6.0d0 6480 else 6481 if (.not.btdb_get(rtdb,'nwpw:wgc_lambda',mt_dbl,1,value)) 6482 > value = 1.0d0 6483 end if 6484 6485 control_wgc_alphabetalambda = value 6486 return 6487 end 6488 6489 6490 6491* *********************************** 6492* * * 6493* * control_mc_step_size * 6494* * * 6495* *********************************** 6496 real*8 function control_mc_step_size() 6497 implicit none 6498 6499#include "bafdecls.fh" 6500#include "btdb.fh" 6501 6502* **** control_rtdb common block **** 6503 integer rtdb 6504 common / control_rtdb1 / rtdb 6505 6506 real*8 value 6507 6508 if (.not.btdb_get(rtdb,'nwpw:mc_step_size',mt_dbl,1,value)) 6509 > value = 0.50d0 6510 6511 control_mc_step_size = value 6512 return 6513 end 6514 6515 6516 6517 6518* *********************************** 6519* * * 6520* * control_mc_volume_step * 6521* * * 6522* *********************************** 6523 real*8 function control_mc_volume_step() 6524 implicit none 6525 6526#include "bafdecls.fh" 6527#include "btdb.fh" 6528 6529* **** control_rtdb common block **** 6530 integer rtdb 6531 common / control_rtdb1 / rtdb 6532 6533 real*8 value 6534 6535 if (.not.btdb_get(rtdb,'nwpw:mc_volume_step',mt_dbl,1,value)) 6536 > value = 0.10d0 6537 6538 control_mc_volume_step = value 6539 return 6540 end 6541 6542 6543 6544 6545* *********************************** 6546* * * 6547* * control_mc_aratio * 6548* * * 6549* *********************************** 6550 real*8 function control_mc_aratio() 6551 implicit none 6552 6553#include "bafdecls.fh" 6554#include "btdb.fh" 6555 6556* **** control_rtdb common block **** 6557 integer rtdb 6558 common / control_rtdb1 / rtdb 6559 6560 real*8 value 6561 6562 if (.not.btdb_get(rtdb,'nwpw:mc_aratio',mt_dbl,1,value)) 6563 > value = 0.234d0 6564 6565 control_mc_aratio = value 6566 return 6567 end 6568 6569 6570* *********************************** 6571* * * 6572* * control_mc_Temperature * 6573* * * 6574* *********************************** 6575 real*8 function control_mc_Temperature() 6576 implicit none 6577 6578#include "bafdecls.fh" 6579#include "btdb.fh" 6580 6581* **** control_rtdb common block **** 6582 integer rtdb 6583 common / control_rtdb1 / rtdb 6584 6585 real*8 value 6586 6587 if (.not.btdb_get(rtdb,'nwpw:mc_temperature',mt_dbl,1,value)) 6588 > value = 298.15d0 6589 6590 control_mc_Temperature = value 6591 return 6592 end 6593 6594 6595* *********************************** 6596* * * 6597* * control_mc_pressure * 6598* * * 6599* *********************************** 6600 real*8 function control_mc_pressure() 6601 implicit none 6602 6603#include "bafdecls.fh" 6604#include "btdb.fh" 6605 6606* **** local variables **** 6607 real*8 autoMbar,autoGPa,autoatm 6608 parameter (autoMbar=294.214239071d0) 6609 parameter (autoGPa=autoMbar*100.0d0) 6610 parameter (autoatm =290.360032539d6) 6611 6612 6613* **** control_rtdb common block **** 6614 integer rtdb 6615 common / control_rtdb1 / rtdb 6616 6617 real*8 value 6618 6619 if (.not.btdb_get(rtdb,'nwpw:mc_pressure',mt_dbl,1,value)) 6620 > value = 1.0d0/autoatm 6621 6622 control_mc_pressure = value 6623 return 6624 end 6625 6626 6627* *********************************** 6628* * * 6629* * control_mc_ddx * 6630* * * 6631* *********************************** 6632 real*8 function control_mc_ddx() 6633 implicit none 6634 6635#include "bafdecls.fh" 6636#include "btdb.fh" 6637 6638* **** control_rtdb common block **** 6639 integer rtdb 6640 common / control_rtdb1 / rtdb 6641 6642 real*8 value 6643 6644 if (.not.btdb_get(rtdb,'nwpw:mc_ddx',mt_dbl,1,value)) 6645 > value = 0.0d0 6646 6647 control_mc_ddx = value 6648 return 6649 end 6650 6651 6652* *********************************** 6653* * * 6654* * control_mc_ddv * 6655* * * 6656* *********************************** 6657 real*8 function control_mc_ddv() 6658 implicit none 6659 6660#include "bafdecls.fh" 6661#include "btdb.fh" 6662 6663* **** control_rtdb common block **** 6664 integer rtdb 6665 common / control_rtdb1 / rtdb 6666 6667 real*8 value 6668 6669 if (.not.btdb_get(rtdb,'nwpw:mc_ddv',mt_dbl,1,value)) 6670 > value = 0.0d0 6671 6672 control_mc_ddv = value 6673 return 6674 end 6675 6676 6677 6678* *********************************** 6679* * * 6680* * control_mc_seed * 6681* * * 6682* *********************************** 6683 integer function control_mc_seed() 6684 implicit none 6685 6686#include "bafdecls.fh" 6687#include "btdb.fh" 6688 6689* **** control_rtdb common block **** 6690 integer rtdb 6691 common / control_rtdb1 / rtdb 6692 6693 integer ivalue 6694 6695 if (.not.btdb_get(rtdb,'nwpw:mc_seed',mt_int,1,ivalue)) 6696 > ivalue = 9484943 6697 6698 control_mc_seed = ivalue 6699 return 6700 end 6701 6702 6703* *********************************** 6704* * * 6705* * control_mc_algorithm * 6706* * * 6707* *********************************** 6708 integer function control_mc_algorithm() 6709 implicit none 6710 6711#include "bafdecls.fh" 6712#include "btdb.fh" 6713 6714* **** control_rtdb common block **** 6715 integer rtdb 6716 common / control_rtdb1 / rtdb 6717 6718 integer ivalue 6719 6720 if (.not.btdb_get(rtdb,'nwpw:mc_algorithm',mt_int,1,ivalue)) 6721 > ivalue = 1 6722 6723 control_mc_algorithm = ivalue 6724 return 6725 end 6726 6727 6728 6729 6730* *********************************** 6731* * * 6732* * control_mc_atom_direction * 6733* * * 6734* *********************************** 6735 subroutine control_mc_atom_direction(mc_atom_direction) 6736 implicit none 6737 real*8 mc_atom_direction(3) 6738 6739#include "bafdecls.fh" 6740#include "btdb.fh" 6741 6742* **** control_rtdb common block **** 6743 integer rtdb 6744 common / control_rtdb1 / rtdb 6745 6746 if (.not.btdb_get(rtdb,'nwpw:mc_atom_direction', 6747 > mt_dbl,3,mc_atom_direction)) then 6748 mc_atom_direction(1) = 1.0d0 6749 mc_atom_direction(2) = 1.0d0 6750 mc_atom_direction(3) = 1.0d0 6751 end if 6752 6753 return 6754 end 6755 6756 6757 6758* *********************************** 6759* * * 6760* * control_mc_ngroups * 6761* * * 6762* *********************************** 6763 subroutine control_mc_ngroups(mc_napply,mc_ngroups,mc_group_size) 6764 implicit none 6765 integer mc_napply,mc_ngroups,mc_group_size 6766 6767#include "bafdecls.fh" 6768#include "btdb.fh" 6769 6770* **** control_rtdb common block **** 6771 integer rtdb 6772 common / control_rtdb1 / rtdb 6773 6774 if (.not.btdb_get(rtdb,'nwpw:mc_napply', 6775 > mt_int,1,mc_napply)) then 6776 mc_napply = 1 6777 end if 6778 if (.not.btdb_get(rtdb,'nwpw:mc_ngroups', 6779 > mt_int,1,mc_ngroups)) then 6780 mc_ngroups = 0 6781 end if 6782 if (.not.btdb_get(rtdb,'nwpw:mc_group_size', 6783 > mt_int,1,mc_group_size)) then 6784 mc_group_size = 0 6785 end if 6786 6787 return 6788 end 6789 6790 6791 6792* *********************************** 6793* * * 6794* * control_mc_groups * 6795* * * 6796* *********************************** 6797 subroutine control_mc_groups(mc_group_start, 6798 > mc_group_end, 6799 > mc_group) 6800 implicit none 6801 integer mc_group_start(*) 6802 integer mc_group_end(*) 6803 integer mc_group(*) 6804 6805#include "bafdecls.fh" 6806#include "btdb.fh" 6807 6808* **** control_rtdb common block **** 6809 integer rtdb 6810 common / control_rtdb1 / rtdb 6811 6812 integer ng,ngs 6813 6814 if (.not.btdb_get(rtdb,'nwpw:mc_ngroups',mt_int,1,ng)) ng=1 6815 if (.not.btdb_get(rtdb,'nwpw:mc_group_size',mt_int,1,ngs)) ngs=1 6816 6817 if (.not.btdb_get(rtdb,'nwpw:mc_group_start', 6818 > mt_int,ng,mc_group_start)) then 6819 call icopy(ng,0,1,mc_group_start,1) 6820 end if 6821 if (.not.btdb_get(rtdb,'nwpw:mc_group_end', 6822 > mt_int,ng,mc_group_end)) then 6823 call icopy(ng,0,1,mc_group_end,1) 6824 end if 6825 if (.not.btdb_get(rtdb,'nwpw:mc_group', 6826 > mt_int,ngs,mc_group)) then 6827 call icopy(ngs,0,1,mc_group,1) 6828 end if 6829 6830 return 6831 end 6832 6833 6834* ************************************ 6835* * * 6836* * control_hess_model * 6837* * * 6838* ************************************ 6839 logical function control_hess_model() 6840 implicit none 6841 6842#include "control.fh" 6843 6844 control_hess_model = hess_model 6845 return 6846 end 6847 6848* ************************************ 6849* * * 6850* * control_hess_filename * 6851* * * 6852* ************************************ 6853 subroutine control_hess_filename(filehess) 6854 implicit none 6855 character*(*) filehess 6856 6857#include "control.fh" 6858#include "bafdecls.fh" 6859#include "btdb.fh" 6860#include "util.fh" 6861 6862* **** control_rtdb common block **** 6863 integer rtdb 6864 common / control_rtdb1 / rtdb 6865 6866 if (.not.btdb_cget(rtdb,'nwpw:hess_model:filename', 6867 > 1,filehess)) then 6868 call util_file_name('hess',.false.,.false.,filehess) 6869 end if 6870 return 6871 end 6872 6873 6874 6875* *********************************** 6876* * * 6877* * control_psp_semicore_small * 6878* * * 6879* *********************************** 6880 logical function control_psp_semicore_small() 6881 implicit none 6882 6883#include "bafdecls.fh" 6884#include "btdb.fh" 6885 6886* **** control_rtdb common block **** 6887 integer rtdb 6888 common / control_rtdb1 / rtdb 6889 6890 logical value 6891 6892 value = .true. 6893 if (.not.btdb_get(rtdb,'nwpw:psp:semicore_small',mt_log,1,value)) 6894 > value = .false. 6895 6896 control_psp_semicore_small = value 6897 return 6898 end 6899 6900* *********************************** 6901* * * 6902* * control_incell * 6903* * * 6904* *********************************** 6905 logical function control_incell() 6906 implicit none 6907 6908#include "bafdecls.fh" 6909#include "btdb.fh" 6910 6911* **** control_rtdb common block **** 6912 integer rtdb 6913 common / control_rtdb1 / rtdb 6914 6915 logical value 6916 6917 value = .true. 6918 if (.not.btdb_get(rtdb,'nwpw:incell',mt_log,1,value)) 6919 > value = .true. 6920 6921 control_incell = value 6922 return 6923 end 6924 6925 6926* *********************************** 6927* * * 6928* * control_fetch_confine * 6929* * * 6930* *********************************** 6931 subroutine control_fetch_confine(symbol,q1,rs,re) 6932 implicit none 6933 character*4 symbol 6934 real*8 q1,rs,re 6935 6936#include "bafdecls.fh" 6937#include "btdb.fh" 6938 6939* **** control_rtdb common block **** 6940 integer rtdb 6941 common / control_rtdb1 / rtdb 6942 6943 logical value 6944 integer l 6945 real*8 rrr(3) 6946 6947 l = index(symbol,' ') - 1 6948 if (.not.btdb_get(rtdb,'nwpw:confine:'//symbol(1:l), 6949 > mt_dbl,3,rrr)) then 6950 6951 rrr(1) = 1.0d0 6952 rrr(2) = 6.0d0 6953 rrr(3) = 14.0d0 6954 end if 6955 q1 = rrr(1) 6956 rs = rrr(2) 6957 re = rrr(3) 6958 return 6959 end 6960 6961 6962* *********************************** 6963* * * 6964* * control_runsocket * 6965* * * 6966* *********************************** 6967 logical function control_runsocket() 6968 implicit none 6969 6970#include "inp.fh" 6971 6972* **** local variables **** 6973 logical value 6974 character*40 stype 6975 6976 value = .false. 6977 call control_socket_type(stype) 6978 if (inp_compare(.false.,'ipi_client',stype)) value = .true. 6979 if (inp_compare(.false.,'unix',stype)) value = .true. 6980 6981 control_runsocket = value 6982 return 6983 end 6984 6985 6986* ************************************ 6987* * * 6988* * control_socket_type * 6989* * * 6990* ************************************ 6991 subroutine control_socket_type(stype) 6992 implicit none 6993 character*(*) stype 6994 6995#include "control.fh" 6996#include "btdb.fh" 6997 6998* **** control_rtdb common block **** 6999 integer rtdb 7000 common / control_rtdb1 / rtdb 7001 7002 if (.not.btdb_cget(rtdb,'nwpw:socket_type', 7003 > 1,stype)) then 7004 stype = "no socket" 7005 end if 7006 return 7007 end 7008 7009* ************************************ 7010* * * 7011* * control_socket_ip * 7012* * * 7013* ************************************ 7014 subroutine control_socket_ip(socket_ip) 7015 implicit none 7016 character*(*) socket_ip 7017 character*40 stype 7018 7019#include "control.fh" 7020#include "btdb.fh" 7021#include "inp.fh" 7022 7023* **** control_rtdb common block **** 7024 integer rtdb 7025 common / control_rtdb1 / rtdb 7026 7027 if (.not.btdb_cget(rtdb,'nwpw:socket_ip', 7028 > 1,socket_ip)) then 7029 call control_socket_type(stype) 7030 if (inp_compare(.false.,'unix',stype)) then 7031 socket_ip = "nwchem" 7032 else 7033 socket_ip = "127.0.0.1:31415" 7034 end if 7035 end if 7036 return 7037 end 7038 7039 7040 7041* ************************************ 7042* * * 7043* * control_vfield_filenames * 7044* * * 7045* ************************************ 7046 logical function control_vfield_filenames(filenames) 7047 implicit none 7048 character*(*) filenames 7049 7050#include "control.fh" 7051#include "btdb.fh" 7052 7053* **** control_rtdb common block **** 7054 integer rtdb 7055 common / control_rtdb1 / rtdb 7056 7057 logical found 7058 7059 found = .false. 7060 if (btdb_cget(rtdb,'nwpw:vfield_filenames',1,filenames)) then 7061 found = .true. 7062 end if 7063 7064 control_vfield_filenames = found 7065 return 7066 end 7067 7068* *********************************** 7069* * * 7070* * control_mult_fixed * 7071* * * 7072* *********************************** 7073 logical function control_mult_fixed() 7074 implicit none 7075 7076#include "bafdecls.fh" 7077#include "btdb.fh" 7078 7079* **** control_rtdb common block **** 7080 integer rtdb 7081 common / control_rtdb1 / rtdb 7082 7083 logical value 7084 7085 value = .true. 7086 if (.not.btdb_get(rtdb,'nwpw:mult_fixed',mt_log,1,value)) 7087 > value = .false. 7088 7089 control_mult_fixed = value 7090 return 7091 end 7092 7093 7094* *********************************** 7095* * * 7096* * control_gas_energy * 7097* * * 7098* *********************************** 7099 real*8 function control_gas_energy() 7100 implicit none 7101 7102#include "bafdecls.fh" 7103#include "btdb.fh" 7104 7105* **** control_rtdb common block **** 7106 integer rtdb 7107 common / control_rtdb1 / rtdb 7108 7109* **** external functions **** 7110 logical nwpw_cosmo_on, nwpw_born_on 7111 external nwpw_cosmo_on, nwpw_born_on 7112 7113 real*8 egas 7114 7115 egas = 0.0d0 7116 if (nwpw_cosmo_on().or.nwpw_born_on()) then 7117 if (.not.btdb_get(rtdb,'nwpw:gas_energy',mt_dbl,1,egas)) 7118 > egas = 0.0d0 7119 end if 7120 7121 control_gas_energy = egas 7122 return 7123 end 7124 7125 7126* *********************************** 7127* * * 7128* * control_gas_energy_set * 7129* * * 7130* *********************************** 7131 logical function control_gas_energy_set(egas) 7132 implicit none 7133 real*8 egas 7134 7135#include "bafdecls.fh" 7136#include "btdb.fh" 7137 7138* **** control_rtdb common block **** 7139 integer rtdb 7140 common / control_rtdb1 / rtdb 7141 7142* **** extertnal functions **** 7143 logical nwpw_cosmo_on,nwpw_born_on 7144 external nwpw_cosmo_on,nwpw_born_on 7145 logical ok 7146 7147 ok=.true. 7148 if (.not.(nwpw_cosmo_on().or.nwpw_born_on())) then 7149 ok=ok.and.btdb_put(rtdb,'nwpw:gas_energy',mt_dbl,1,egas) 7150 end if 7151 control_gas_energy_set=ok 7152 return 7153 end 7154 7155 7156* *********************************** 7157* * * 7158* * control_H1_it_in * 7159* * * 7160* *********************************** 7161 integer function control_H1_it_in() 7162 implicit none 7163 7164#include "bafdecls.fh" 7165#include "btdb.fh" 7166 7167* **** control_rtdb common block **** 7168 integer rtdb 7169 common / control_rtdb1 / rtdb 7170 7171 integer it_in 7172 7173 if (.not.btdb_get(rtdb,'nwpw:H1_it_in',mt_int,1,it_in)) 7174 > it_in = 50 7175 7176 control_H1_it_in = it_in 7177 return 7178 end 7179 7180* *********************************** 7181* * * 7182* * control_H1_it_out * 7183* * * 7184* *********************************** 7185 integer function control_H1_it_out() 7186 implicit none 7187 7188#include "bafdecls.fh" 7189#include "btdb.fh" 7190 7191* **** control_rtdb common block **** 7192 integer rtdb 7193 common / control_rtdb1 / rtdb 7194 7195 integer it_out 7196 7197 if (.not.btdb_get(rtdb,'nwpw:H1_it_out',mt_int,1,it_out)) 7198 > it_out = 1 7199 7200 control_H1_it_out = it_out 7201 return 7202 end 7203 7204 7205* *********************************** 7206* * * 7207* * control_H1_it_ortho * 7208* * * 7209* *********************************** 7210 integer function control_H1_it_ortho() 7211 implicit none 7212 7213#include "bafdecls.fh" 7214#include "btdb.fh" 7215 7216* **** control_rtdb common block **** 7217 integer rtdb 7218 common / control_rtdb1 / rtdb 7219 7220 integer it_ortho 7221 7222 if (.not.btdb_get(rtdb,'nwpw:H1_it_ortho',mt_int,1,it_ortho)) 7223 > it_ortho = -1 7224 7225 control_H1_it_ortho = it_ortho 7226 return 7227 end 7228 7229* *********************************** 7230* * * 7231* * control_CI_maxit_sweeps * 7232* * * 7233* *********************************** 7234 integer function control_CI_maxit_sweeps() 7235 implicit none 7236 7237#include "bafdecls.fh" 7238#include "btdb.fh" 7239 7240* **** control_rtdb common block **** 7241 integer rtdb 7242 common / control_rtdb1 / rtdb 7243 7244 integer it_out 7245 7246 if (.not.btdb_get(rtdb,'nwpw:CI_maxit_sweeps',mt_int,1,it_out)) 7247 > it_out = 3 7248 7249 control_CI_maxit_sweeps = it_out 7250 return 7251 end 7252 7253 7254* *********************************** 7255* * * 7256* * control_CI_maxit_orb * 7257* * * 7258* *********************************** 7259 integer function control_CI_maxit_orb() 7260 implicit none 7261 7262#include "bafdecls.fh" 7263#include "btdb.fh" 7264 7265* **** control_rtdb common block **** 7266 integer rtdb 7267 common / control_rtdb1 / rtdb 7268 7269 integer it_out 7270 7271 if (.not.btdb_get(rtdb,'nwpw:CI_maxit_orb',mt_int,1,it_out)) 7272 > it_out = 120 7273 7274 control_CI_maxit_orb = it_out 7275 return 7276 end 7277 7278* *********************************** 7279* * * 7280* * control_CI_maxit_linesearch * 7281* * * 7282* *********************************** 7283 integer function control_CI_maxit_linesearch() 7284 implicit none 7285 7286#include "bafdecls.fh" 7287#include "btdb.fh" 7288 7289* **** control_rtdb common block **** 7290 integer rtdb 7291 common / control_rtdb1 / rtdb 7292 7293 integer it_out 7294 7295 if (.not.btdb_get(rtdb,'nwpw:CI_maxit_linesearch', 7296 > mt_int,1,it_out)) 7297 > it_out = 15 7298 7299 control_CI_maxit_linesearch = it_out 7300 return 7301 end 7302 7303 7304* *********************************** 7305* * * 7306* * control_use_director * 7307* * * 7308* *********************************** 7309 logical function control_use_director() 7310 implicit none 7311 7312#include "bafdecls.fh" 7313#include "btdb.fh" 7314 7315* **** control_rtdb common block **** 7316 integer rtdb 7317 common / control_rtdb1 / rtdb 7318 7319 logical value 7320 7321 value = .true. 7322 if (.not.btdb_get(rtdb,'nwpw:use_director',mt_log,1,value)) 7323 > value = .false. 7324 7325 control_use_director = value 7326 return 7327 end 7328 7329 7330* *********************************** 7331* * * 7332* * control_unset_excited_ne * 7333* * * 7334* *********************************** 7335 subroutine control_unset_excited_ne() 7336 implicit none 7337 7338#include "bafdecls.fh" 7339#include "btdb.fh" 7340 7341* **** control_rtdb common block **** 7342 integer rtdb 7343 common / control_rtdb1 / rtdb 7344 7345 logical value 7346 integer ne(2) 7347 7348 ne(1) = 0 7349 ne(2) = 0 7350 value = btdb_put(rtdb,'nwpw:excited_ne',mt_int,2,ne) 7351 7352 !value = rtdb_delete(rtdb,'nwpw:excited_ne') 7353 7354 return 7355 end 7356 7357 7358* *********************************** 7359* * * 7360* * control_Efield_on * 7361* * * 7362* *********************************** 7363 logical function control_Efield_on() 7364 implicit none 7365 7366#include "bafdecls.fh" 7367#include "btdb.fh" 7368 7369* **** control_rtdb common block **** 7370 integer rtdb 7371 common / control_rtdb1 / rtdb 7372 7373 logical value 7374 7375 value = .true. 7376 if (.not.btdb_get(rtdb,'nwpw:efield',mt_log,1,value)) 7377 > value = .false. 7378 7379 control_Efield_on = value 7380 return 7381 end 7382 7383 7384* ********************************************* 7385* * * 7386* * control_single_precision_on * 7387* * * 7388* ********************************************* 7389 logical function control_single_precision_on() 7390 implicit none 7391 7392#include "control.fh" 7393 7394 control_single_precision_on = single_precision_on 7395 return 7396 end 7397