TABLE OF CONTENTS
ABINIT/m_multibinit_main [ Modules ]
NAME
m_multibinit_main
FUNCTION
Main routine MULTIBINIT.
COPYRIGHT
Copyright (C) 1999-2024 ABINIT group (AM, hexu) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
(main routine)
OUTPUT
(main routine)
NOTES
Should be 1 moved to somewhere else 2 be replaced with the new implementation multibinit_main2.
SOURCE
28 #if defined HAVE_CONFIG_H 29 #include "config.h" 30 #endif 31 32 #include "abi_common.h" 33 34 ! FIXME: This module is temporarily here. It should be removed once we have the new lattice mover. 35 ! Move the main calculation part into a subroutine 36 ! TODO: And this subroutine should be replaced by a one liner then: call multibinit_manager%run_all(). 37 ! TODO: The module is here because it uses lv95_drive mover_effpot (very hacky). 38 ! TODO: The mover_effpot is in 95_drive due to that it use some >78 level module 39 ! TODO: which should be changed when we implement the new lattice mover. 40 module m_multibinit_driver 41 use defs_basis 42 use defs_abitypes 43 use m_xmpi 44 use m_xomp 45 use m_abicore 46 use m_errors 47 48 use m_effective_potential 49 use m_fit_polynomial_coeff 50 use m_opt_effpot 51 use m_multibinit_dataset 52 use m_effective_potential_file 53 use m_scup_dataset 54 !use m_spin_model, only: spin_model_t 55 use m_abihist 56 57 use m_build_info, only: abinit_version 58 use m_multibinit_manager, only: mb_manager_t 59 use m_multibinit_main2, only: multibinit_main2 60 61 use m_mover_effpot, only : mover_effpot 62 #if defined DEV_MS_SCALEUP 63 use scup_global 64 #endif 65 !use m_generate_training_set, only : generate_training_set 66 use m_compute_anharmonics, only : compute_anharmonics 67 use m_init10, only : init10, postfix_fnames 68 use m_parser, only : instrng 69 use m_fstrings, only : replace, inupper 70 use m_dtset, only : chkvars 71 implicit none
m_multbinit_main/multibinit_main [ Functions ]
NAME
multibinit_main
FUNCTION
The main function of multibinit
INPUTS
filnam: The filenames from the files file. 17 files in total.
OUTPUT
SOURCE
88 subroutine multibinit_main(input_path, filnam, dry_run) 89 character(len=fnlen), intent(inout) :: input_path 90 character(len=fnlen), intent(inout) :: filnam(18) 91 integer, intent(in) :: dry_run 92 type(multibinit_dtset_type), target :: inp 93 type(effective_potential_type) :: reference_effective_potential 94 type(abihist) :: hist, hist_tes 95 96 !type(spin_model_t) :: spin_model 97 character(len=strlen) :: string, raw_string 98 character(len=500) :: message 99 character(len=fnlen) :: name 100 character(len=fnlen) :: sys_fname 101 102 integer :: filetype,ii,lenstr,iiter,niter 103 integer :: natom,nph1l,nrpt,ntypat 104 integer :: option 105 logical :: need_analyze_anh_pot,need_prt_files 106 107 ! Whether the "new" MULTIBNIT framework should be used. 108 logical :: need_new_multibinit 109 ! MS 110 ! temporary variables for testing SCALE-UP with Multibinit 111 !Variable to pass to effpot_evaluate routine of multibinit 112 !To declare evaluation of electronice model 113 logical :: elec_eval 114 #if defined DEV_MS_SCALEUP 115 !Variables needed to call SCALE-UP 116 logical :: err_init_elec 117 logical*1 :: needlattice = .FALSE. 118 logical*1 :: needelectrons = .TRUE. 119 logical*1 :: didi = .FALSE. 120 logical*1 :: harm_der = .FALSE. 121 logical*1 :: initorbocc = .FALSE. 122 logical*1 :: ismagnetic = .FALSE. 123 logical*1 :: istddft = .FALSE. 124 logical*4 :: printgeom = .FALSE. 125 logical*4 :: printeigv = .FALSE. 126 logical*4 :: printeltic = .FALSE. 127 logical*4 :: printorbocc = .FALSE. 128 integer :: ksamp(3) 129 real*8 :: tcharge 130 #endif 131 !TEST_AM 132 ! integer :: natom_sp 133 ! real(dp),allocatable :: dynmat(:,:,:,:,:) 134 !TEST_AM 135 !****************************************************************** 136 137 integer :: master, my_rank, comm, nproc, ierr 138 logical :: iam_master 139 140 141 !MPI variables 142 master = 0 143 comm = xmpi_world 144 nproc = xmpi_comm_size(comm) 145 my_rank = xmpi_comm_rank(comm) 146 iam_master = (my_rank == master) 147 148 !Read the input file, only to the the name of the file which contains the ddb file or xml file 149 ! for getting the number of atoms. 150 option=1 151 if (iam_master) then 152 call instrng (filnam(1),lenstr,option,strlen,string,raw_string) 153 !To make case-insensitive, map characters to upper case: 154 call inupper(string(1:lenstr)) 155 !Check whether the string only contains valid keywords 156 call chkvars(string) 157 end if 158 call xmpi_bcast(string,master, comm, ierr) 159 call xmpi_bcast(raw_string,master, comm, ierr) 160 call xmpi_bcast(lenstr,master, comm, ierr) 161 !To automate a maximum calculation, multibinit reads the number of atoms 162 !in the file (ddb or xml). If DDB file is present in input, the ifc calculation 163 !will be initilaze array to the maximum of atoms (natifc=natom,atifc=1,natom...) in invars10 164 165 166 !Read the input file assuming natom=1 so that the invars10 can work. 167 168 !call invars10(inp,lenstr,natom,string) 169 if(trim(filnam(3)) /='') then 170 sys_fname=filnam(3) 171 else 172 call invars_multibinit_filenames(string, lenstr, sys_fname=sys_fname) 173 end if 174 175 176 177 ! read the reference structure to get natom 178 if (iam_master) then 179 write(message, '(6a)' )' Read the information in the reference structure in ',ch10,& 180 & '-',trim(sys_fname),ch10,' to initialize the multibinit input' 181 call wrtout(ab_out,message,'COLL') 182 call wrtout(std_out,message,'COLL') 183 end if 184 185 call effective_potential_file_getDimSystem(sys_fname,comm,natom,ntypat,nph1l,nrpt) 186 !call effective_potential_file_getDimSystem(filnam(3),natom,ntypat,nph1l,nrpt) 187 188 189 ! read the input again to use the right natom 190 call invars10(inp,lenstr,natom,string) 191 call postfix_fnames(input_path, filnam, inp) 192 193 need_new_multibinit= inp%spin_dynamics > 0 .or. inp%lwf_dynamics > 0 .or. inp%dynamics >= 100 194 195 if (iam_master) then 196 if(need_new_multibinit) then 197 ABI_ERROR("The new MULTINIT mode should be enabled with --F03 option. ") 198 end if 199 ! Echo the inputs to console and main output file 200 call outvars_multibinit(inp,std_out) 201 call outvars_multibinit(inp,ab_out) 202 end if 203 204 if(dry_run/=0) then 205 call wrtout([std_out, ab_out], "Multibinit in dry_run mode. Exiting after input parser") 206 call xmpi_end() 207 !goto 100 208 endif 209 210 ! Read the model (from DDB or XML) 211 call effective_potential_file_read(filnam(3),reference_effective_potential,inp,comm) 212 213 214 if(filnam(4)/=''.and.filnam(4)/='no') then 215 call effective_potential_file_getType(filnam(4),filetype) 216 if(filetype==3.or.filetype==23) then 217 call effective_potential_file_read(filnam(4),reference_effective_potential,inp,comm) 218 else 219 write(message,'(a,(80a),3a)') ch10,('=',ii=1,80),ch10,ch10,& 220 & ' There is no specific file for the coefficients from polynomial fitting' 221 call wrtout(ab_out,message,'COLL') 222 call wrtout(std_out,message,'COLL') 223 end if 224 else 225 if(inp%ncoeff/=0) then 226 write(message, '(5a)' )& 227 & 'ncoeff is specified in the input but,',ch10,& 228 & 'there is no file for the coefficients ',ch10,& 229 & 'Action: add coefficients.xml file' 230 ABI_ERROR(message) 231 232 else 233 write(message,'(a,(80a),3a)') ch10,('=',ii=1,80),ch10,ch10,& 234 & ' There is no file for the coefficients from polynomial fitting' 235 call wrtout(ab_out,message,'COLL') 236 call wrtout(std_out,message,'COLL') 237 end if 238 end if 239 240 !**************************************************************************************** 241 !SCALE UP Initialize the electronic model (If scale-up is available) 242 !**************************************************************************************** 243 elec_eval = .FALSE. 244 245 #if defined DEV_MS_SCALEUP 246 if(inp%scup_dtset%scup_elec_model)then 247 write(message,'(a,(80a),4a)') ch10,('=',ii=1,80),ch10,ch10,& 248 ' Initializing Electronic Model with SCALE-UP',ch10 249 call wrtout(ab_out,message,'COLL') 250 call wrtout(std_out,message,'COLL') 251 252 !Set Variables 253 elec_eval = .TRUE. 254 ksamp = inp%scup_dtset%scup_ksamp 255 tcharge = inp%scup_dtset%scup_tcharge 256 if(inp%scup_dtset%scup_ismagnetic)ismagnetic=.TRUE. 257 if(inp%scup_dtset%scup_istddft)istddft=.TRUE. 258 if(inp%scup_dtset%scup_initorbocc)initorbocc=.TRUE. 259 260 ! Call to Scale-Up 261 err_init_elec = global_init_model(filnam(3),inp%ncell,needlattice,needelectrons,didi,& 262 & harm_der,tcharge,ksamp,ismagnetic,istddft,initorbocc) 263 264 !Set Print variables 265 if(inp%scup_dtset%scup_printgeom)printgeom=.TRUE. 266 if(inp%scup_dtset%scup_printeigv)printeigv=.TRUE. 267 if(inp%scup_dtset%scup_printeltic)printeltic=.TRUE. 268 if(inp%scup_dtset%scup_printorbocc)printorbocc=.TRUE. 269 270 !Set Print Parameters within scaleup 271 call global_set_print_parameters(printgeom,printeigv,printeltic,printorbocc,& 272 & inp%scup_dtset%scup_printbands) 273 274 !Set SCF controling variables (values contain defaults, if not specified in the input) 275 call global_set_scf_parameters(inp%scup_dtset%scup_scfmixing,inp%scup_dtset%scup_scfthresh,& 276 & inp%scup_dtset%scup_smearing,inp%scup_dtset%scup_maxscfstep,& 277 & inp%scup_dtset%scup_startpulay,inp%scup_dtset%scup_freezden) 278 279 280 !Create kpath if printbands=true and pass it to SCALE UP 281 if(inp%scup_dtset%scup_printbands)then 282 283 call scup_kpath_new(inp%scup_dtset%scup_speck,& 284 & reference_effective_potential%supercell%rprimd,& 285 & inp%scup_dtset%scup_ndivsm,inp%scup_dtset%scup_kpath) 286 call scup_kpath_print(inp%scup_dtset%scup_kpath) 287 288 call global_set_print_bands(inp%scup_dtset%scup_printbands,& 289 & inp%scup_dtset%scup_nspeck,inp%scup_dtset%scup_kpath%ndivs,& 290 & inp%scup_dtset%scup_speck) 291 endif 292 endif 293 #endif 294 295 !**************************************************************************************** 296 ! Compute the third order derivative with finite differences 297 !**************************************************************************************** 298 if (inp%strcpling > 0) then 299 call compute_anharmonics(reference_effective_potential,filnam,inp,comm) 300 end if 301 !**************************************************************************************** 302 303 ! If needed, fit the anharmonic part and compute the confinement potential 304 !**************************************************************************************** 305 if (inp%fit_coeff/=0.or.inp%confinement==2.or.inp%bound_model/=0 .or. inp%opt_effpot/=0) then 306 307 if(iam_master) then 308 ! Read the MD file 309 write(message,'(a,(80a),7a)')ch10,('=',ii=1,80),ch10,ch10,& 310 & '-Reading the training-set file :',ch10,& 311 & '-',trim(filnam(5)),ch10 312 313 call wrtout(std_out,message,'COLL') 314 call wrtout(ab_out,message,'COLL') 315 if(filnam(5)/=''.and.filnam(5)/='no')then 316 call effective_potential_file_readMDfile(filnam(5),hist,option=inp%ts_option) 317 if (hist%mxhist == 0)then 318 write(message, '(5a)' )& 319 & 'The trainig-set ',trim(filnam(5)),' file is not correct ',ch10,& 320 & 'Action: add training-set file' 321 ABI_ERROR(message) 322 end if 323 else 324 if (inp%fit_coeff/=0) then 325 write(message, '(3a)' )& 326 & 'There is no training-set file to fit the lattice model ',ch10,& 327 & 'Action: add trainings-set-file' 328 ABI_ERROR(message) 329 else if (inp%bound_model/=0) then 330 write(message, '(3a)' )& 331 & 'There is no training-set file to bound the model ',ch10,& 332 & 'Action: add training-set file ' 333 ABI_ERROR(message) 334 else if(inp%confinement==2) then 335 write(message, '(3a)' )& 336 & 'There is no training-set file to compute the confinement',ch10,& 337 & 'Action: add training-set file ' 338 ABI_ERROR(message) 339 else if(inp%opt_effpot==2) then 340 write(message, '(3a)' )& 341 & 'There is no training-set file to optimize the latice model',ch10,& 342 & 'Action: add training-set file ' 343 ABI_ERROR(message) 344 end if 345 end if 346 end if 347 ! MPI BROADCAST the history of the MD 348 call abihist_bcast(hist,master,comm) 349 ! Map the hist in order to be consistent with the supercell into reference_effective_potential 350 call effective_potential_file_mapHistToRef(reference_effective_potential,hist,comm) 351 end if 352 353 !TEST_AM 354 ! call effective_potential_checkDEV(reference_effective_potential,hist,size(hist%xred,2),hist%mxhist) 355 ! stop 356 !TEST_AM 357 358 !Generate the confinement polynome (not working yet) 359 if(inp%confinement/=0)then 360 option=inp%confinement 361 select case(option) 362 case(1) 363 call effective_potential_setConfinement(inp%conf_cutoff_disp,inp%conf_cutoff_strain,& 364 & reference_effective_potential,inp%conf_power_fact_disp,& 365 & inp%conf_power_fact_strain,inp%conf_power_disp,& 366 & inp%conf_power_disp,inp%conf_power_strain,& 367 & need_confinement=.TRUE.) 368 369 write(message,'(a,(80a),3a)') ch10,('=',ii=1,80),ch10,ch10,& 370 & ' The confinement potential is active.' 371 call wrtout(ab_out,message,'COLL') 372 call wrtout(std_out,message,'COLL') 373 374 case(2) 375 write(message,'(a,(80a),3a)') ch10,('=',ii=1,80),ch10,ch10,& 376 & ' The confinement potential is computed from the MD file and actived.' 377 call wrtout(ab_out,message,'COLL') 378 call wrtout(std_out,message,'COLL') 379 380 end select 381 end if 382 383 384 !Fit the coeff 385 if (inp%fit_coeff/=0)then 386 option=inp%fit_coeff 387 if(hist%mxhist >0)then 388 if (option==-1)then 389 ! option == -1 390 ! Print the file in the specific format for the script of carlos 391 ! Born_Charges 392 ! Dielectric_Tensor 393 ! harmonic.xml 394 ! Reference_structure 395 ! Strain_Tensor 396 ! symmetry_operations (only cubic) 397 if (iam_master) then 398 call fit_polynomial_printSystemFiles(reference_effective_potential,hist) 399 end if 400 else if (option==1.or.option==2)then 401 ! option = 1 402 if(inp%fit_iatom/=0)then 403 call fit_polynomial_coeff_fit(reference_effective_potential,& 404 & inp%fit_bancoeff,inp%fit_fixcoeff,hist,inp%fit_generateCoeff,& 405 & inp%fit_rangePower,inp%fit_nbancoeff,inp%fit_ncoeff,& 406 & inp%fit_nfixcoeff,inp%fit_nimposecoeff,inp%fit_imposecoeff,& 407 & option,comm,cutoff_in=inp%fit_cutoff,& 408 & max_power_strain=inp%fit_SPC_maxS,initialize_data=inp%fit_initializeData==1,& 409 & fit_tolMSDF=inp%fit_tolMSDF,fit_tolMSDS=inp%fit_tolMSDS,fit_tolMSDE=inp%fit_tolMSDE,& 410 & fit_tolMSDFS=inp%fit_tolMSDFS,fit_tolGF=inp%fit_tolGF,& 411 & verbose=.true.,positive=.false.,& 412 & anharmstr=inp%fit_anhaStrain==1,& 413 & spcoupling=inp%fit_SPCoupling==1,prt_anh=inp%analyze_anh_pot,& 414 & fit_iatom=inp%fit_iatom,prt_files=.TRUE.,fit_on=inp%fit_on,sel_on=inp%sel_on,& 415 & fit_factors=inp%fit_factors,prt_GF_csv=inp%prt_GF_csv,dispterms=inp%fit_dispterms==1) 416 else 417 if (inp%fit_ncoeff_per_iatom/=0)then 418 if (mod(inp%fit_ncoeff,inp%fit_ncoeff_per_iatom) /= 0)then 419 write(message,'(2a,I3,2a,I3,3a)') ch10,& 420 & 'fit_ncoeff_per_iatom = ', inp%fit_ncoeff_per_iatom,ch10,& 421 & 'is not a divider of fit_ncoeff = ', inp%fit_ncoeff,ch10,& 422 & 'Action: Change fit_ncoeff and/or fit_ncoeff_per_iatom',ch10 423 ABI_ERROR(message) 424 endif 425 niter = inp%fit_ncoeff/inp%fit_ncoeff_per_iatom 426 if (mod(niter,reference_effective_potential%crystal%nirredat) /= 0)then 427 write(message,'(2a,I3,2a,I3,2a,I3,3a)') ch10,& 428 & 'fit_ncoeff_per_iatom = ', inp%fit_ncoeff_per_iatom,ch10,& 429 & 'times the number of irreducible atoms = ',reference_effective_potential%crystal%nirredat,ch10,& 430 & 'is not a divider of fit_ncoeff = ', inp%fit_ncoeff,ch10,& 431 & 'Action: Change fit_ncoeff and/or fit_ncoeff_per_iatom',ch10 432 ABI_ERROR(message) 433 endif 434 niter = niter/reference_effective_potential%crystal%nirredat 435 else if (inp%fit_ncoeff_per_iatom == 0)then 436 if (mod(inp%fit_ncoeff,reference_effective_potential%crystal%nirredat) /= 0)then 437 write(message,'(2a,I3,2a,I3,3a)') ch10,& 438 & 'The number of irreducible atoms = ',reference_effective_potential%crystal%nirredat,ch10,& 439 & 'is not a divider of fit_ncoeff = ', inp%fit_ncoeff,ch10,& 440 & 'Action: Change fit_ncoeff',ch10 441 ABI_ERROR(message) 442 endif 443 inp%fit_ncoeff_per_iatom = inp%fit_ncoeff/reference_effective_potential%crystal%nirredat 444 niter = 1 445 endif 446 write(message,'(a,(80a),7a,I3,3a,I3,3a,I3,3a,I3,2a)') ch10,('=',ii=1,80),ch10,ch10,& 447 & ' Starting Fit Iterations ',ch10,& 448 & ' ----------------------- ',ch10,& 449 & ' Select in total fit_ncoeff = ', inp%fit_ncoeff,' coefficients',ch10,& 450 & ' In ', niter,' iterations',ch10,& 451 & ' Over ', reference_effective_potential%crystal%nirredat, ' irreducible atoms',ch10,& 452 & ' Selecting ', inp%fit_ncoeff_per_iatom, ' coefficients per atom in each iteration',ch10 453 call wrtout(std_out,message,'COLL') 454 call wrtout(ab_out,message,'COLL') 455 need_prt_files=.FALSE. 456 do iiter=1,niter 457 write(message,'(a,(80a),3a,I3,a,I3,2a)') ch10,('-',ii=1,80),ch10,ch10,& 458 & ' Start Iteration (',iiter,'/',niter,')',ch10 459 call wrtout(std_out,message,'COLL') 460 call wrtout(ab_out,message,'COLL') 461 do ii=1,reference_effective_potential%crystal%nirredat 462 if(ii == reference_effective_potential%crystal%nirredat .and. iiter==niter)need_prt_files=.TRUE. 463 if(ii > 1 .or. iiter > 1)inp%fit_nfixcoeff = -1 464 call fit_polynomial_coeff_fit(reference_effective_potential,& 465 & inp%fit_bancoeff,inp%fit_fixcoeff,hist,inp%fit_generateCoeff,& 466 & inp%fit_rangePower,inp%fit_nbancoeff,inp%fit_ncoeff_per_iatom,& 467 & inp%fit_nfixcoeff,inp%fit_nimposecoeff,inp%fit_imposecoeff,& 468 & option,comm,cutoff_in=inp%fit_cutoff,& 469 & max_power_strain=inp%fit_SPC_maxS,initialize_data=inp%fit_initializeData==1,& 470 & fit_tolMSDF=inp%fit_tolMSDF,fit_tolMSDS=inp%fit_tolMSDS,fit_tolMSDE=inp%fit_tolMSDE,& 471 & fit_tolMSDFS=inp%fit_tolMSDFS,fit_tolGF=inp%fit_tolGF,& 472 & verbose=.true.,positive=.false.,& 473 & anharmstr=inp%fit_anhaStrain==1,& 474 & spcoupling=inp%fit_SPCoupling==1,prt_anh=inp%analyze_anh_pot,& 475 & fit_iatom=reference_effective_potential%crystal%irredatindx(ii),& 476 & prt_files=need_prt_files,fit_on=inp%fit_on,sel_on=inp%sel_on,& 477 & fit_factors=inp%fit_factors,prt_GF_csv=inp%prt_GF_csv,dispterms=inp%fit_dispterms==1) 478 enddo 479 enddo 480 endif 481 end if 482 else 483 write(message, '(3a)' )& 484 & 'There is no step in the MD file ',ch10,& 485 & 'Action: add correct MD file' 486 ABI_ERROR(message) 487 end if 488 end if 489 490 491 !try to bound the model with mover_effpot 492 !we need to use the molecular dynamics 493 if(inp%bound_model>0.and.inp%bound_model<=2)then 494 call mover_effpot(inp,filnam,reference_effective_potential,-1*inp%bound_model,comm,hist=hist) 495 !Marcus: New option for bound_model: use optimize routine for generting specific high order terms 496 elseif(inp%bound_model == 3)then 497 write(message,'(a,(80a),4a)')ch10,('=',ii=1,80),ch10,ch10,& 498 & 'Bound Process 3: Generate equivalent high order terms',ch10 499 call wrtout(std_out,message,'COLL') 500 call wrtout(ab_out,message,'COLL') 501 502 call opt_effpotbound(reference_effective_potential,inp%bound_rangePower,hist, inp%bound_EFS,& 503 & inp%bound_factors,inp%bound_penalty,comm) 504 505 end if 506 507 !**************************************************************************************** 508 ! OPTIMIZE SECTION, Optimize selected coefficients of effective potential while 509 ! keeping the others constant 510 !**************************************************************************************** 511 512 if(inp%opt_effpot == 1)then 513 write(message,'(a,(80a),4a)')ch10,('=',ii=1,80),ch10,ch10,& 514 & 'Optimizing Effective Potential',ch10 515 516 call wrtout(std_out,message,'COLL') 517 call wrtout(ab_out,message,'COLL') 518 519 need_analyze_anh_pot = .FALSE. 520 if(inp%analyze_anh_pot == 1) need_analyze_anh_pot = .TRUE. 521 522 call opt_effpot(reference_effective_potential,inp%opt_ncoeff,inp%opt_coeff,hist,inp%opt_on,& 523 & inp%opt_factors,comm,print_anh=need_analyze_anh_pot) 524 end if 525 526 527 528 !**************************************************************************************** 529 ! TEST SECTION test effective potential with regard to test-set 530 !**************************************************************************************** 531 if(inp%test_effpot == 1)then 532 if(iam_master) then 533 ! Read the test-set .nc file 534 write(message,'(a,(80a),9a)')ch10,('=',ii=1,80),ch10,ch10,& 535 & 'TEST - SET Option',ch10,& 536 & '-Reading the test-set file :',ch10,& 537 & '-',trim(filnam(6)),ch10 538 539 call wrtout(std_out,message,'COLL') 540 call wrtout(ab_out,message,'COLL') 541 if(filnam(6)/=''.and.filnam(6)/='no')then 542 call effective_potential_file_readMDfile(filnam(6),hist_tes,option=inp%ts_option) 543 if (hist_tes%mxhist == 0)then 544 write(message, '(5a)' )& 545 & 'The test-set ',trim(filnam(6)),' file is empty ',ch10,& 546 & 'Action: add non-empty test-set' 547 ABI_ERROR(message) 548 end if 549 else 550 write(message, '(3a)' )& 551 & 'There is no test-set file ',ch10,& 552 & 'Action: add test-set file' 553 ABI_ERROR(message) 554 end if 555 end if 556 ! MPI BROADCAST the history of the MD 557 call abihist_bcast(hist_tes,master,comm) 558 ! Map the hist in order to be consistent with the supercell into reference_effective_potential 559 call effective_potential_file_mapHistToRef(reference_effective_potential,hist_tes,comm) 560 ! Initialize if to print anharmonic contribution to energy or not 561 need_analyze_anh_pot = .FALSE. 562 if(inp%analyze_anh_pot == 1) need_analyze_anh_pot = .TRUE. 563 ! Call to test routine 564 call fit_polynomial_coeff_testEffPot(reference_effective_potential,hist_tes,master,comm,& 565 & print_anharmonic=need_analyze_anh_pot,scup_dtset=inp%scup_dtset,& 566 & prt_ph=inp%test_prt_ph) 567 568 569 570 end if ! End if(inp%test_effpot == 1)then 571 572 !TEST_AM 573 !Effective Hamiltonian, compute the energy for given patern 574 ! call mover_effpot(inp,filnam,reference_effective_potential,-2,comm,hist=hist) 575 !TEST_AM 576 577 !**************************************************************************************** 578 579 !**************************************************************************************** 580 !Print the effective potential system + coefficients (only master CPU) 581 if(iam_master) then 582 if (inp%prt_model >= 1) then 583 write(message, '(a,(80a),a)' ) ch10,& 584 & ('=',ii=1,80) 585 call wrtout(ab_out,message,'COLL') 586 call wrtout(std_out,message,'COLL') 587 !name = replace(trim(filnam(2)),".out","") 588 ! Assume new .abo convention 589 name = replace(trim(filnam(2)),".abo","") 590 call effective_potential_writeXML(reference_effective_potential,inp%prt_model,filename=name,& 591 & prt_dipdip=inp%dipdip_prt==1) 592 else if (inp%prt_model == -2)then 593 ! NetCDF case, in progress 594 name = trim(filnam(2))//"_sys.nc" 595 call effective_potential_writeNETCDF(reference_effective_potential,1,filename=name) 596 end if 597 end if 598 !**************************************************************************************** 599 600 !TEST_AM SECTION 601 ! Print the Phonon dos/spectrum 602 ! if(inp%prt_phfrq > 0) then 603 ! call effective_potential_printPDOS(reference_effective_potential,filnam(2),& 604 !& inp%ncell,inp%nph1l,inp%prt_phfrq,inp%qph1l) 605 ! end if 606 607 !Intialisation of the effective potential type 608 ! call effective_potential_file_read(filnam(4),reference_effective_potential,inp,comm) 609 ! name = "test.xml" 610 ! call effective_potential_writeXML(reference_effective_potential,1,filename=name) 611 ! just for TEST 612 ! if(inp%prt_phfrq > 0) then 613 ! natom_sp = reference_effective_potential%supercell%natom_supercell 614 ! ABI_MALLOC(dynmat,(2,3,natom_sp,3,natom_sp)) 615 ! call effective_potential_effpot2dynmat(dynmat,inp%delta_df,reference_effective_potential,& 616 ! & reference_effective_potential%supercell%natom_supercell,& 617 ! & int(reference_effective_potential%supercell%qphon),3) 618 619 ! ABI_FREE(dynmat) 620 ! end if 621 ! end if 622 !TEST_AM SECTION 623 624 625 ! Run lattice dynamics (relaxation or molecular dynamics, most of abinits ionmovs are allowed) 626 !**************************************************************************************** 627 if(inp%dynamics>=1) then 628 call mover_effpot(inp,filnam,reference_effective_potential,inp%dynamics,comm) 629 end if 630 631 !**************************************************************************************** 632 633 634 635 !Free the effective_potential and dataset 636 !**************************************************************************************** 637 call effective_potential_free(reference_effective_potential) 638 call multibinit_dtset_free(inp) 639 call abihist_free(hist) 640 call abihist_free(hist_tes) 641 !**************************************************************************************** 642 643 end subroutine multibinit_main