TABLE OF CONTENTS
ABINIT/mover_effpot [ Functions ]
NAME
mover_effpot
FUNCTION
this routine is driver for using mover with effective potential
COPYRIGHT
Copyright (C) 1998-2017 ABINIT group (AM) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
inp = input of multibinit effective_potential = effective potential of the reference structure option = flag for the option: -1 = Bound the anharmonic part 12 = NVT simulation 13 = NPT simulation
OUTPUT
SIDE EFFECTS
NOTES
PARENTS
multibinit
CHILDREN
alloc_copy,destroy_mpi_enreg,destroy_results_gs,dtset_free effective_potential_setcoeffs,effective_potential_setsupercell fit_polynomial_coeff_fit,fit_polynomial_coeff_getpositive,generelist init_results_gs,mover,polynomial_coeff_free,polynomial_coeff_getnorder polynomial_coeff_init,polynomial_coeff_setcoefficient polynomial_coeff_writexml,scfcv_destroy,wrtout,xcart2xred,xmpi_barrier xred2xcart
SOURCE
44 #if defined HAVE_CONFIG_H 45 #include "config.h" 46 #endif 47 48 #include "abi_common.h" 49 50 subroutine mover_effpot(inp,filnam,effective_potential,option,comm,hist) 51 52 use defs_basis 53 use defs_abitypes 54 use m_profiling_abi 55 use defs_datatypes 56 use m_errors 57 use m_abimover 58 use m_build_info 59 use m_scf_history 60 use defs_wvltypes 61 use m_xmpi 62 use m_abimover 63 use m_phonons 64 use m_strain 65 use m_effective_potential_file 66 use m_supercell 67 use m_psps 68 use m_args_gs 69 70 use m_geometry, only : xcart2xred, xred2xcart 71 use m_multibinit_dataset, only : multibinit_dtset_type 72 use m_effective_potential,only : effective_potential_type 73 use m_fit_polynomial_coeff, only : polynomial_coeff_writeXML 74 use m_fit_polynomial_coeff, only : fit_polynomial_coeff_fit,genereList 75 use m_fit_polynomial_coeff, only : fit_polynomial_coeff_getPositive,fit_polynomial_coeff_getCoeffBound 76 use m_electronpositron, only : electronpositron_type 77 use m_polynomial_coeff,only : polynomial_coeff_getNorder 78 ! use m_pawang, only : pawang_type, pawang_free 79 ! use m_pawrad, only : pawrad_type, pawrad_free 80 ! use m_pawtab, only : pawtab_type, pawtab_nullify, pawtab_free 81 ! use m_pawxmlps, only : paw_setup, ipsp2xml, rdpawpsxml, & 82 !& paw_setup_copy, paw_setup_free, getecutfromxml 83 use m_dtset, only : dtset_free 84 use m_abihist, only : abihist 85 use m_ifc 86 use m_ewald 87 use m_mpinfo, only : init_mpi_enreg,destroy_mpi_enreg 88 use m_copy , only : alloc_copy 89 use m_electronpositron, only : electronpositron_type 90 use m_scfcv, only : scfcv_t, scfcv_run,scfcv_destroy 91 use m_results_gs, only : results_gs_type,init_results_gs,destroy_results_gs 92 93 !This section has been created automatically by the script Abilint (TD). 94 !Do not modify the following lines by hand. 95 #undef ABI_FUNC 96 #define ABI_FUNC 'mover_effpot' 97 use interfaces_14_hidewrite 98 use interfaces_95_drive, except_this_one => mover_effpot 99 !End of the abilint section 100 101 implicit none 102 103 !Arguments -------------------------------- 104 !scalar 105 integer, intent(in) :: option,comm 106 !array 107 type(multibinit_dtset_type),intent(in) :: inp 108 type(effective_potential_type),intent(inout) :: effective_potential 109 character(len=fnlen),intent(in) :: filnam(15) 110 type(abihist),optional,intent(inout):: hist 111 !Local variables------------------------------- 112 !scalar 113 integer :: icoeff_bound,ii 114 !integer :: iexit,initialized 115 integer :: jj,kk,nproc,ncoeff,nmodels,ncoeff_bound,ncoeff_bound_tot,ncoeff_max 116 integer :: model_bound,model_ncoeffbound,my_rank 117 !integer :: mtypalch,,npsp,paw_size,type 118 !integer,save :: paw_size_old=-1 119 real(dp):: cutoff,freq_q,freq_b,qmass,bmass 120 logical :: iam_master 121 integer, parameter:: master=0 122 logical :: verbose,writeHIST 123 !real(dp) :: cpui 124 !character(len=6) :: codvsn 125 126 !TEST_AM 127 ! integer :: ia,mu,rand_seed = 5 128 ! real(dp):: mass_ia,rescale_vel,sum_mass,v2gauss 129 !TEST_AM 130 ! Set array dimensions 131 character(len=500) :: message 132 type(MPI_type),target :: mpi_enreg 133 type(dataset_type),target :: dtset 134 type(scfcv_t) :: scfcv_args 135 type(datafiles_type),target :: dtfil 136 integer,target :: zero_integer 137 type(ab_xfh_type) :: ab_xfh 138 type(results_gs_type),target :: results_gs 139 type(pseudopotential_type),target :: psps 140 !arrays 141 !no_abirules 142 integer :: sc_size(3),sc_size_TS(3) 143 integer,pointer :: indsym(:,:,:) 144 integer,allocatable :: listcoeff(:),listcoeff_bound(:,:),list_tmp(:),list_bound(:,:) 145 integer,allocatable :: isPositive(:) 146 integer,allocatable :: symrel(:,:,:) 147 !integer,allocatable :: npwtot(:) 148 real(dp) :: acell(3) 149 !real(dp) :: ecut_tmp(3,2,10) 150 real(dp),allocatable :: amass(:) ,coeff_values(:,:) 151 real(dp),pointer :: rhog(:,:),rhor(:,:) 152 real(dp),allocatable :: tnons(:,:) 153 real(dp),allocatable :: xred(:,:),xred_old(:,:),xcart(:,:) 154 real(dp),allocatable :: fred(:,:),fcart(:,:) 155 real(dp),allocatable :: vel(:,:) 156 real(dp) :: vel_cell(3,3),rprimd(3,3) 157 type(polynomial_coeff_type),dimension(:),allocatable :: coeffs_all,coeffs_tmp,coeffs_bound 158 character(len=fnlen) :: filename 159 !character(len=fnlen) :: filename_psp(3) 160 type(electronpositron_type),pointer :: electronpositron 161 ! type(pspheader_type),allocatable :: pspheads(:) 162 ! type(pawrad_type),allocatable :: pawrad(:) 163 ! type(pawtab_type),allocatable :: pawtab(:) 164 ! type(args_gs_type) :: args_gs 165 !type(wvl_data) :: wvl 166 !type(pawang_type) :: pawang 167 !type(scf_history_type) :: scf_history 168 169 !****************************************************************** 170 171 !MPI variables 172 nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 173 iam_master = (my_rank == master) 174 175 write(message, '(a,(80a),a)') ch10,& 176 & ('=',ii=1,80),ch10 177 call wrtout(ab_out,message,'COLL') 178 call wrtout(std_out,message,'COLL') 179 180 !******************************************************************* 181 ! 1 Generate and check supercell for the dynamics 182 !******************************************************************* 183 184 !a new supercell is compute 185 !Initialisaton of variable 186 if(option == -1.or.option == -2) then 187 ! Bound process option 188 sc_size(:) = inp%fit_boundCell 189 else if(option == -3) then 190 ! Heff option 191 sc_size(:) = (/1,1,1/) 192 else 193 ! Normal dynamics 194 sc_size(:) = inp%n_cell 195 end if 196 197 if(option/=0)then 198 199 acell = one 200 rprimd = effective_potential%crystal%rprimd 201 202 ABI_ALLOCATE(xred,(3,effective_potential%crystal%natom)) 203 ABI_ALLOCATE(xcart,(3,effective_potential%crystal%natom)) 204 205 ! convert new xcart 206 call xcart2xred(effective_potential%crystal%natom,effective_potential%crystal%rprimd,& 207 & effective_potential%crystal%xcart,xred) 208 call xred2xcart(effective_potential%crystal%natom, rprimd, xcart, xred) 209 ! Generate supercell for the simulation 210 call effective_potential_setSupercell(effective_potential,comm,n_cell=sc_size) 211 212 ABI_DEALLOCATE(xred) 213 ABI_DEALLOCATE(xcart) 214 215 !*************************************************************** 216 !1 Convert some parameters into the structures used by mover.F90 217 !*************************************************************** 218 219 ! Free dtset 220 call dtset_free(dtset) 221 222 ! Set mpi_eng 223 mpi_enreg%comm_cell = comm 224 mpi_enreg%me = my_rank 225 226 ! Set the abinit dataset for mover with fake values 227 ! Scalar 228 dtset%dmft_entropy = 0 229 dtset%nctime = inp%nctime ! NetCdf TIME between output of molecular dynamics informations 230 dtset%delayperm = 0 ! DELAY between trials to PERMUTE atoms 231 dtset%dilatmx = 1.0 ! DILATation : MaXimal value 232 dtset%diismemory = 8 ! Direct Inversion in the Iterative Subspace MEMORY 233 dtset%friction = 0.0001 ! internal FRICTION coefficient 234 dtset%goprecon = 0 ! Geometry Optimization PREconditioner equations 235 dtset%istatr = 0 ! Integer for STATus file SHiFT 236 dtset%jellslab = 0 ! include a JELLium SLAB in the cell 237 dtset%mqgrid = 0 ! Maximum number of Q-space GRID points for pseudopotentials 238 dtset%mqgriddg = 0 ! Maximum number of Q-wavevectors for the 1-dimensional GRID 239 ! for the Double Grid in PAW 240 dtset%mdwall = 10000 ! Molecular Dynamics WALL location 241 dtset%ntypalch = 0 ! Number of TYPe of atoms that are "ALCHemical" 242 dtset%natom = effective_potential%supercell%natom 243 dtset%ntypat = effective_potential%crystal%ntypat 244 dtset%npspalch = effective_potential%crystal%ntypat 245 dtset%nconeq = 0 ! Number of CONstraint EQuations 246 dtset%noseinert = 1.d-5 ! NOSE INERTia factor 247 dtset%nnos = inp%nnos ! Number of nose masses Characteristic 248 dtset%nsym = 1 ! Number of SYMmetry operations 249 dtset%prtxml = 0 ! print the xml 250 dtset%signperm = 1 ! SIGN of PERMutation potential 251 dtset%strprecon = 1 ! STRess PRECONditioner 252 dtset%tolmxf = 2.0d-5 253 dtset%tsmear = 0.009500446 ! 254 dtset%vis = 100 ! VIScosity 255 dtset%usewvl = 0 ! 256 dtset%useylm = 0 ! 257 258 ! if(option == -3) then 259 ! write(message,'(a)')' Read the DDB file to fill the dtset array' 260 ! call wrtout(std_out,message,"COLL") 261 ! ! Copy real informtions from the ddb 262 ! call effective_potential_file_getType(filnam(3),type) 263 ! if(type /= 1) then 264 ! write(message, '(5a)' )& 265 ! & ' You need to provide DDB file in the input to compute ahnarmonic',ch10,& 266 ! & ' part of effective Hamiltionian',ch10,& 267 ! & 'Action: add DDB file in the inputs' 268 ! MSG_BUG(message) 269 ! end if 270 ! call ddb_to_dtset(comm, dtset,filnam(3),psps) 271 ! ABI_ALLOCATE(dtset%kptns,(3,dtset%nkpt)) 272 ! dtset%kptns(:,:) = dtset%kpt(:,:) 273 ! ABI_ALLOCATE(dtset%istwfk,(dtset%nkpt)) 274 ! dtset%istwfk(:) = 1 275 276 ! else 277 ! Need to init some values 278 ABI_ALLOCATE(symrel,(3,3,dtset%nsym)) 279 symrel = 1 280 call alloc_copy(symrel,dtset%symrel) 281 ABI_ALLOCATE(tnons,(3,dtset%nsym)) 282 tnons = zero 283 call alloc_copy(tnons,dtset%tnons) 284 call alloc_copy(effective_potential%supercell%typat,dtset%typat) 285 call alloc_copy(effective_potential%crystal%znucl,dtset%znucl) 286 ABI_DEALLOCATE(symrel) 287 ABI_DEALLOCATE(tnons) 288 ! end if 289 290 !array 291 ABI_ALLOCATE(dtset%iatfix,(3,dtset%natom)) ! Indices of AToms that are FIXed 292 dtset%iatfix = 0 293 dtset%goprecprm(:) = zero !Geometry Optimization PREconditioner PaRaMeters equations 294 ABI_ALLOCATE(dtset%prtatlist,(dtset%natom)) !PRinT by ATom LIST of ATom 295 dtset%prtatlist(:) = 0 296 ABI_ALLOCATE(dtset%mixalch_orig,(dtset%npspalch,dtset%ntypalch,1)) 297 dtset%mixalch_orig(:,:,:)=zero 298 if(option > 0)then 299 verbose = .TRUE. 300 writeHIST = .TRUE. 301 dtset%dtion = inp%dtion ! Delta Time for IONs 302 dtset%ionmov = inp%dynamics ! Number for the dynamic 303 dtset%ntime = inp%ntime ! Number of TIME steps 304 dtset%optcell = inp%optcell ! OPTimize the CELL shape and dimensions Characteristic 305 dtset%restartxf = inp%restartxf ! RESTART from (X,F) history 306 dtset%mdtemp(1) = inp%temperature !Molecular Dynamics Temperatures 307 dtset%mdtemp(2) = inp%temperature !Molecular Dynamics Temperatures 308 dtset%strtarget(1:6) = -1 * inp%strtarget(1:6) / 29421.033d0 ! STRess TARGET 309 else if(option == -1.or.option == -2) then 310 ! Set default for the fit 311 verbose = .false. 312 writeHIST = .false. 313 dtset%restartxf = 0 ! RESTART from (X,F) history 314 dtset%dtion = 100 ! Delta Time for IONs 315 dtset%ionmov = 13 ! Number for the dynamic 316 dtset%ntime = inp%fit_boundStep ! Number of TIME steps 317 dtset%optcell = 2 ! OPTimize the CELL shape and dimensions Characteristic 318 dtset%mdtemp(1) = inp%fit_boundTemp !Molecular Dynamics Temperatures 319 dtset%mdtemp(2) = inp%fit_boundTemp !Molecular Dynamics Temperatures 320 dtset%strtarget(1:6) = zero 321 end if 322 323 ! Set the barostat and thermonstat if ionmov == 13 324 if(dtset%ionmov == 13)then 325 326 ! Select frequency of the barostat as a function of temperature 327 ! For small temperature, we need huge barostat and inversely 328 ! if(dtset%mdtemp(1) <= 10) then 329 ! freq_q = 0.002 330 ! freq_b = 0.0002 331 ! else if(dtset%mdtemp(1) <= 50) then 332 ! freq_q = 0.02 333 ! freq_b = 0.002 334 ! else if(dtset%mdtemp(1) > 50.and.dtset%mdtemp(1) < 300) then 335 ! freq_q = 0.1 336 ! freq_b = 0.01 337 ! else 338 ! freq_q = 0.2 339 ! freq_b = 0.02 340 ! end if 341 342 !TEST_AM 343 freq_q = 0.1 344 freq_b = 0.01 345 !TEST_AM 346 347 qmass = dtset%natom* kb_THzK * dtset%mdtemp(1) / (freq_q**2) 348 bmass = dtset%natom* kb_THzK * dtset%mdtemp(1) / (freq_b**2) 349 350 if(dtset%nnos==0) then 351 dtset%nnos = 1 352 ABI_ALLOCATE(dtset%qmass,(dtset%nnos)) 353 dtset%qmass(:) = qmass 354 write(message,'(3a,F20.1,a)')& 355 & ' WARNING: nnos is set to zero in the input',ch10,& 356 & ' value by default for qmass: ',dtset%qmass(:),ch10 357 if(verbose)call wrtout(std_out,message,"COLL") 358 else 359 ABI_ALLOCATE(dtset%qmass,(dtset%nnos)) ! Q thermostat mass 360 dtset%qmass(:) = inp%qmass(:) 361 end if 362 if (abs(inp%bmass) < tol10) then 363 dtset%bmass = bmass 364 write(message,'(3a,F20.4,a)')& 365 & ' WARNING: bmass is set to zero in the input',ch10,& 366 & ' value by default for bmass: ',dtset%bmass,ch10 367 if(verbose)call wrtout(std_out,message,"COLL") 368 else 369 dtset%bmass = inp%bmass ! Barostat mass 370 end if 371 end if 372 373 374 ! set psps 375 psps%useylm = dtset%useylm 376 ! if(option == -3)then 377 ! mtypalch = 0 378 ! npsp = dtset%ntypat 379 ! call psps_free(psps) 380 ! filename_psp(1) = "/home/alex/calcul/psp/Sr.LDA_PW-JTH.xml" 381 ! filename_psp(2) = "/home/alex/calcul/psp/Ti.LDA_PW-JTH.xml" 382 ! filename_psp(3) = "/home/alex/calcul/psp/O.LDA_PW-JTH.xml" 383 ! ABI_DATATYPE_ALLOCATE(pspheads,(npsp)) 384 ! call inpspheads(filename_psp,npsp,pspheads,ecut_tmp) 385 ! call psps_init_global(mtypalch, npsp, psps, pspheads) 386 ! call psps_init_from_dtset(dtset, 1, psps, pspheads) 387 ! end if 388 389 ! ! The correct dimension of pawrad/tab is ntypat. In case of alchemical psps 390 ! ! pawrad/tab(ipsp) is invoked with ipsp<=npsp. So, in order to avoid any problem, 391 ! ! declare pawrad/tab at paw_size=max(ntypat,npsp). 392 ! paw_size=0;if (psps%usepaw==1) paw_size=max(dtset%ntypat,dtset%npsp) 393 ! if (paw_size/=paw_size_old) then 394 ! if (paw_size_old/=-1) then 395 ! call pawrad_free(pawrad) 396 ! call pawtab_free(pawtab) 397 ! ABI_DATATYPE_DEALLOCATE(pawrad) 398 ! ABI_DATATYPE_DEALLOCATE(pawtab) 399 ! end if 400 ! ABI_DATATYPE_ALLOCATE(pawrad,(paw_size)) 401 ! ABI_DATATYPE_ALLOCATE(pawtab,(paw_size)) 402 ! call pawtab_nullify(pawtab) 403 ! paw_size_old=paw_size 404 ! end if 405 406 ! set args_gs 407 ! if (option == -3 then 408 ! call args_gs_init(args_gs, & 409 ! & effective_potential%crystal%amu(:),dtset%mixalch_orig(:,:,1),& 410 ! & dtset%dmatpawu(:,:,:,:,1),dtset%upawu(:,1),dtset%jpawu(:,1),& 411 ! & dtset%rprimd_orig(:,:,1)) 412 ! ABI_ALLOCATE(npwtot,(dtset%nkpt)) 413 ! end if 414 ! initialisation of results_gs 415 call init_results_gs(dtset%natom,1,results_gs) 416 417 ! Set the pointers of scfcv_args 418 zero_integer = 0 419 scfcv_args%dtset => dtset 420 ABI_ALLOCATE(indsym,(4,dtset%nsym,dtset%natom)) 421 indsym = 0 422 scfcv_args%indsym => indsym 423 scfcv_args%mpi_enreg => mpi_enreg 424 scfcv_args%ndtpawuj => zero_integer 425 scfcv_args%results_gs => results_gs 426 scfcv_args%psps => psps 427 ! Set other arguments of the mover.F90 routines 428 429 ABI_ALLOCATE(amass,(dtset%natom)) 430 ! Assign masses to each atom (for MD) 431 do jj = 1,dtset%natom 432 amass(jj)=amu_emass*& 433 & effective_potential%crystal%amu(effective_potential%supercell%typat(jj)) 434 end do 435 ! Set the dffil structure 436 dtfil%filnam_ds(1:2)=filnam(1:2) 437 dtfil%filnam_ds(3)="" 438 dtfil%filnam_ds(4)=filnam(2) 439 dtfil%filstat='_STATUS' 440 nullify (electronpositron) 441 ABI_ALLOCATE(rhog,(2,1)) 442 ABI_ALLOCATE(rhor,(2,1)) 443 444 ! Initialize xf history (should be put in inwffil) 445 ! Not yet implemented for ionmov 2 3 10 11 22 (memory problem...) 446 ! ab_xfh%mxfh=(ab_xfh%nxfh-dtset%restartxf+1)+dtset%ntime+5 447 ab_xfh%nxfh = 0 448 ab_xfh%mxfh = 1 449 ABI_ALLOCATE(ab_xfh%xfhist,(3,dtset%natom+4,2,ab_xfh%mxfh)) 450 if (any((/2,3,10,11,22/)==dtset%ionmov)) then 451 write(message, '(3a)' )& 452 & ' This dynamics can not be used with effective potential',ch10,& 453 & 'Action: correct dynamics input' 454 MSG_BUG(message) 455 end if 456 457 !*************************************************************** 458 !2 initialization of the structure for the dynamics 459 !*************************************************************** 460 461 if (allocated(dtset%rprimd_orig)) then 462 ABI_DEALLOCATE(dtset%rprimd_orig) 463 end if 464 ABI_ALLOCATE(dtset%rprimd_orig,(3,3,1)) 465 dtset%rprimd_orig(:,:,1) = effective_potential%supercell%rprimd 466 467 acell(1) = dtset%rprimd_orig(1,1,1) 468 acell(2) = dtset%rprimd_orig(2,2,1) 469 acell(3) = dtset%rprimd_orig(3,3,1) 470 471 ABI_ALLOCATE(xred,(3,dtset%natom)) 472 ABI_ALLOCATE(xred_old,(3,dtset%natom)) 473 ABI_ALLOCATE(vel,(3,dtset%natom)) 474 ABI_ALLOCATE(fred,(3,dtset%natom)) 475 ABI_ALLOCATE(fcart,(3,dtset%natom)) 476 477 call xcart2xred(dtset%natom,effective_potential%supercell%rprimd,& 478 & effective_potential%supercell%xcart,xred) 479 480 xred_old = xred 481 vel_cell(:,:) = zero 482 vel(:,:) = zero 483 484 !********************************************************* 485 !4 Call main routine for the bound process, 486 ! monte carlo / molecular dynamics / project 487 !********************************************************* 488 if(option > 0)then 489 !************************************************************* 490 ! call mover in case of NPT or NVT simulation 491 !************************************************************* 492 write(message, '((80a),3a)' ) ('-',ii=1,80), ch10,& 493 & '-Monte Carlo / Molecular Dynamics ',ch10 494 call wrtout(ab_out,message,'COLL') 495 call wrtout(std_out,message,'COLL') 496 call mover(scfcv_args,ab_xfh,acell,amass,dtfil,electronpositron,& 497 & rhog,rhor,dtset%rprimd_orig,vel,vel_cell,xred,xred_old,& 498 & effective_potential=effective_potential,verbose=verbose,writeHIST=writeHIST) 499 500 else if(option== -1.or.option==-2)then 501 !************************************************************* 502 ! Try to bound the model 503 !************************************************************* 504 write(message, '((80a),4a)' ) ('-',ii=1,80), ch10,& 505 & ' Try to bound the model',ch10,' Check if the model is bounded or not' 506 call wrtout(ab_out,message,'COLL') 507 call wrtout(std_out,message,'COLL') 508 509 ! Try the model 510 call mover(scfcv_args,ab_xfh,acell,amass,dtfil,electronpositron,& 511 & rhog,rhor,dtset%rprimd_orig,vel,vel_cell,xred,xred_old,& 512 & effective_potential=effective_potential,verbose=verbose,writeHIST=writeHIST) 513 514 write(message, '(a)' ) ' => The model' 515 if(effective_potential%anharmonics_terms%bounded)then 516 write(message, '(2a)' ) trim(message),' is bound' 517 call wrtout(std_out,message,'COLL') 518 call wrtout(ab_out,message,'COLL') 519 else 520 write(message, '(2a)' ) trim(message),' is not bound' 521 call wrtout(std_out,message,'COLL') 522 523 524 if(option==-2)then 525 526 ! Fill the list for the fixcoeff input of the fit_polynomial_coeff_fit routine 527 ! Store the number of coefficients before adding other coeff for the bounding 528 ncoeff = effective_potential%anharmonics_terms%ncoeff 529 ABI_ALLOCATE(listcoeff,(ncoeff)) 530 do ii=1,ncoeff 531 listcoeff(ii)=ii 532 end do 533 534 write(message, '(2a)')ch10,' Generate the list of addionnal terms with the fit process...' 535 call wrtout(std_out,message,'COLL') 536 537 ! Get the additional coeff 538 call fit_polynomial_coeff_fit(effective_potential,(/0/),listcoeff,hist,1,& 539 & inp%fit_boundPower,0,inp%fit_boundTerm,ncoeff,1,comm,cutoff_in=inp%fit_boundCutoff,& 540 & max_power_strain=2,verbose=.true.,positive=.true.,spcoupling=inp%fit_SPCoupling==1,& 541 & anharmstr=inp%fit_anhaStrain==1,only_even_power=.true.) 542 543 ! Store the max number of coefficients after the fit process 544 ncoeff_max = effective_potential%anharmonics_terms%ncoeff 545 ! Store all the coefficients in coeffs_all 546 ABI_DATATYPE_ALLOCATE(coeffs_all,(ncoeff_max)) 547 ABI_DATATYPE_ALLOCATE(coeffs_tmp,(ncoeff_max)) 548 do ii=1,ncoeff_max 549 call polynomial_coeff_init(& 550 & effective_potential%anharmonics_terms%coefficients(ii)%coefficient,& 551 & effective_potential%anharmonics_terms%coefficients(ii)%nterm,& 552 & coeffs_all(ii),& 553 & effective_potential%anharmonics_terms%coefficients(ii)%terms,& 554 & effective_potential%anharmonics_terms%coefficients(ii)%name,& 555 & check=.false.) 556 end do 557 558 model_ncoeffbound = 0 559 do ii=1,ncoeff_max-ncoeff 560 561 write(message, '(5a,I0,a)')ch10,'--',ch10,' Try to bound the model ',& 562 & 'with ', ii,' additional term' 563 if(ii>1)write(message,'(2a)') trim(message),'s' 564 call wrtout(std_out,message,'COLL') 565 566 ! Copy the new model in coeffs_tmp(jj) 567 ! Free the coeffs_tmp array before 568 do jj=1,ncoeff_max 569 call polynomial_coeff_free(coeffs_tmp(jj)) 570 end do 571 572 do jj=1,ncoeff+ii 573 call polynomial_coeff_init(& 574 & coeffs_all(jj)%coefficient,& 575 & coeffs_all(jj)%nterm,& 576 & coeffs_tmp(jj),& 577 & coeffs_all(jj)%terms,& 578 & coeffs_all(jj)%name,& 579 & check=.false.) 580 end do 581 582 model_ncoeffbound = ii 583 584 ! Reset the simulation and set the coefficients of the model 585 call effective_potential_setCoeffs(coeffs_tmp(1:ncoeff+ii),effective_potential,ncoeff+ii) 586 call fit_polynomial_coeff_fit(effective_potential,(/0/),(/0/),hist,0,(/0,0/),0,0,& 587 & -1,1,comm,verbose=.true.,positive=.false.) 588 call effective_potential_setSupercell(effective_potential,comm,n_cell=sc_size) 589 dtset%rprimd_orig(:,:,1) = effective_potential%supercell%rprimd 590 acell(1) = dtset%rprimd_orig(1,1,1) 591 acell(2) = dtset%rprimd_orig(2,2,1) 592 acell(3) = dtset%rprimd_orig(3,3,1) 593 call xcart2xred(dtset%natom,effective_potential%supercell%rprimd,& 594 & effective_potential%supercell%xcart,xred) 595 xred_old = xred 596 vel_cell(:,:) = zero 597 vel(:,:) = zero 598 fred(:,:) = zero 599 fcart(:,:) = zero 600 601 ! Run mover to check if the model is bound 602 call mover(scfcv_args,ab_xfh,acell,amass,dtfil,electronpositron,& 603 & rhog,rhor,dtset%rprimd_orig,vel,vel_cell,xred,xred_old,& 604 & effective_potential=effective_potential,verbose=verbose,writeHIST=writeHIST) 605 if(.not.effective_potential%anharmonics_terms%bounded)then 606 write(message, '(2a)' ) ' => The model is not bounded' 607 else 608 write(message, '(2a)' ) ' => The model is bounded' 609 end if 610 call wrtout(std_out,message,'COLL') 611 612 ! Exit if the model is bounded 613 if(effective_potential%anharmonics_terms%bounded) exit 614 615 end do 616 617 else 618 ! Get the list of possible coefficients to bound the model 619 cutoff = zero 620 do ii=1,3 621 cutoff = cutoff + effective_potential%crystal%rprimd(ii,ii) 622 end do 623 cutoff = cutoff / 3.0 624 625 ! call fit_polynomial_coeff_getCoeffBound(effective_potential,coeffs_bound,& 626 !& hist,ncoeff_bound,comm,verbose=.true.) 627 628 !TEST_AM! 629 sc_size_TS = (/2,2,2/) 630 call polynomial_coeff_getNorder(coeffs_bound,effective_potential%crystal,cutoff,& 631 & ncoeff_bound,ncoeff_bound_tot,inp%fit_boundPower,inp%fit_boundPower(2),2,sc_size_TS,& 632 & comm,anharmstr=inp%fit_anhaStrain==1,& 633 & spcoupling=inp%fit_SPCoupling==1,verbose=.false.,distributed=.false.,& 634 & only_even_power=.true.,only_odd_power=.false.) 635 636 if(iam_master)then 637 filename=trim(filnam(2))//"_boundcoeff.xml" 638 call polynomial_coeff_writeXML(coeffs_bound,ncoeff_bound,filename=filename,newfile=.true.) 639 end if 640 ! wait 641 call xmpi_barrier(comm) 642 ! Store all the initial coefficients 643 ncoeff = effective_potential%anharmonics_terms%ncoeff 644 ABI_DATATYPE_ALLOCATE(coeffs_all,(ncoeff+ncoeff_bound)) 645 do ii=1,ncoeff 646 call polynomial_coeff_init(effective_potential%anharmonics_terms%coefficients(ii)%coefficient,& 647 & effective_potential%anharmonics_terms%coefficients(ii)%nterm,& 648 & coeffs_all(ii),& 649 & effective_potential%anharmonics_terms%coefficients(ii)%terms,& 650 & effective_potential%anharmonics_terms%coefficients(ii)%name,& 651 & check=.false.) 652 end do 653 654 do ii=1,ncoeff_bound 655 call polynomial_coeff_init(coeffs_bound(ii)%coefficient,& 656 & coeffs_bound(ii)%nterm,& 657 & coeffs_all(ncoeff+ii),& 658 & coeffs_bound(ii)%terms,& 659 & coeffs_bound(ii)%name,& 660 & check=.false.) 661 end do 662 663 ! Copy the fixed coefficients from the model (without bound coeff) 664 ncoeff = effective_potential%anharmonics_terms%ncoeff 665 ABI_DATATYPE_ALLOCATE(coeffs_tmp,(ncoeff+ncoeff_bound)) 666 do ii=1,ncoeff 667 call polynomial_coeff_init(effective_potential%anharmonics_terms%coefficients(ii)%coefficient,& 668 & effective_potential%anharmonics_terms%coefficients(ii)%nterm,& 669 & coeffs_tmp(ii),& 670 & effective_potential%anharmonics_terms%coefficients(ii)%terms,& 671 & effective_potential%anharmonics_terms%coefficients(ii)%name,& 672 & check=.false.) 673 end do 674 675 ncoeff_max = ncoeff+ncoeff_bound 676 ABI_ALLOCATE(listcoeff,(ncoeff_max)) 677 listcoeff = 0 678 do jj=1,ncoeff 679 listcoeff(jj) = jj 680 end do 681 682 model_bound = 0 683 model_ncoeffbound = 0 684 685 do ii=2,inp%fit_boundTerm 686 ! Compute the number of possible combination 687 nmodels = 1 688 ABI_ALLOCATE(list_bound,(nmodels,ii)) 689 ABI_ALLOCATE(list_tmp,(ii)) 690 list_bound = 0; list_tmp = 0; kk = 0; jj = 1 691 692 ! Generate the list of possible combinaison 1st count 693 call genereList(kk,jj,ii,ncoeff_bound,list_tmp,list_bound,nmodels,.false.) 694 nmodels = kk 695 696 write(message, '(5a,I0,a,I0,a)')ch10,'--',ch10,' Try to bound the model ',& 697 & 'with ', ii,' additional positive terms (',nmodels,') possibilities' 698 call wrtout(std_out,message,'COLL') 699 700 ! allocate and generate combinaisons 701 ABI_DEALLOCATE(list_bound) 702 ABI_DEALLOCATE(list_tmp) 703 ABI_ALLOCATE(coeff_values,(nmodels,ncoeff+ii)) 704 ABI_ALLOCATE(listcoeff_bound,(nmodels,ncoeff+ii)) 705 ABI_ALLOCATE(list_bound,(nmodels,ii)) 706 ABI_ALLOCATE(list_tmp,(ii)) 707 ABI_ALLOCATE(isPositive,(nmodels)) 708 list_bound = 0; listcoeff_bound = 0; list_tmp = 0; isPositive = 0; kk = 0; jj = 1 709 call genereList(kk,jj,ii,ncoeff_bound,list_tmp,list_bound,nmodels,.true.) 710 ! Generate the models 711 do jj=1,nmodels 712 listcoeff_bound(jj,1:ncoeff) = listcoeff(1:ncoeff) 713 listcoeff_bound(jj,ncoeff+1:ncoeff+ii) = list_bound(jj,:) + ncoeff 714 end do 715 716 ! Reset the simulation 717 call effective_potential_setCoeffs(coeffs_all,effective_potential,ncoeff+ncoeff_bound) 718 call fit_polynomial_coeff_getPositive(effective_potential,hist,coeff_values,& 719 & isPositive,listcoeff_bound,ncoeff+ii,ncoeff,nmodels,comm,verbose=.false.) 720 if(all(isPositive == 0)) then 721 write(message, '(5a,I0,a)')ch10,'--',ch10,' No possible model ',& 722 & 'with ', ii,' additional terms found' 723 call wrtout(std_out,message,'COLL') 724 else 725 726 do jj=1,nmodels 727 if(isPositive(jj) == 1 .and. all(abs(coeff_values(jj,:)) < 1.0E5)) then 728 write(message, '(2a,I0,a)') ch10,' The model number ',jj,' [' 729 do kk=1,ncoeff+ii 730 if(kk<ncoeff+ii)then 731 write(message, '(a,I0,a)') trim(message),listcoeff_bound(jj,kk),',' 732 else 733 write(message, '(a,I0)') trim(message),listcoeff_bound(jj,kk) 734 end if 735 end do 736 write(message, '(2a)') trim(message),'] is positive' 737 call wrtout(std_out,message,'COLL') 738 write(message, '(2a,I0,a)') ' Check if the model ',& 739 & 'number ', jj,' is bounded...' 740 call wrtout(std_out,message,'COLL') 741 742 ! Set the coefficients of the model 743 do kk=1,ncoeff+ii 744 if(kk<=ncoeff)then 745 ! just set the values of the coefficient 746 call polynomial_coeff_setCoefficient(coeff_values(jj,kk),coeffs_tmp(kk)) 747 write(message, '(a,I0,a,ES19.10,2a)') ' Set the value of the coefficient ',kk,& 748 & ' =>',coeff_values(jj,kk),' ',trim(coeffs_tmp(kk)%name) 749 call wrtout(std_out,message,'COLL') 750 751 else 752 ! Set the good coefficient 753 icoeff_bound = listcoeff_bound(jj,kk)-ncoeff ! need to remove ncoeff value 754 write(message, '(a,I0,a,I0,a,ES19.10,2a)')& 755 & ' Set the value of the coefficient ',kk,' (',icoeff_bound,& 756 & ') =>',coeff_values(jj,kk),& 757 & ' ',trim(coeffs_bound(icoeff_bound)%name) 758 call wrtout(std_out,message,'COLL') 759 call polynomial_coeff_free(coeffs_tmp(kk)) 760 call polynomial_coeff_init(coeff_values(jj,kk),& 761 & coeffs_bound(icoeff_bound)%nterm,& 762 & coeffs_tmp(kk),& 763 & coeffs_bound(icoeff_bound)%terms,& 764 & coeffs_bound(icoeff_bound)%name,& 765 & check=.false.) 766 767 end if 768 end do 769 770 ! Reset the simulation and set the coefficients of the model 771 call effective_potential_setCoeffs(coeffs_tmp(1:ncoeff+ii),effective_potential,& 772 & ncoeff+ii) 773 call fit_polynomial_coeff_fit(effective_potential,(/0/),(/0/),hist,0,(/0,0/),1,0,& 774 & -1,1,comm,verbose=.false.,positive=.false.) 775 call effective_potential_setSupercell(effective_potential,comm,n_cell=sc_size) 776 dtset%rprimd_orig(:,:,1) = effective_potential%supercell%rprimd 777 acell(1) = dtset%rprimd_orig(1,1,1) 778 acell(2) = dtset%rprimd_orig(2,2,1) 779 acell(3) = dtset%rprimd_orig(3,3,1) 780 call xcart2xred(dtset%natom,effective_potential%supercell%rprimd,& 781 & effective_potential%supercell%xcart,xred) 782 xred_old = xred 783 vel_cell(:,:) = zero 784 vel(:,:) = zero 785 fred(:,:) = zero 786 fcart(:,:) = zero 787 788 ! Run mover 789 call mover(scfcv_args,ab_xfh,acell,amass,dtfil,electronpositron,& 790 & rhog,rhor,dtset%rprimd_orig,vel,vel_cell,xred,xred_old,& 791 & effective_potential=effective_potential,verbose=verbose,writeHIST=writeHIST) 792 793 if(.not.effective_potential%anharmonics_terms%bounded)then 794 write(message, '(2a)' ) ' => The model is not bounded' 795 else 796 write(message, '(2a)' ) ' => The model is bounded' 797 end if 798 call wrtout(std_out,message,'COLL') 799 ! Exit if the model is bounded 800 if(effective_potential%anharmonics_terms%bounded) then 801 model_bound = jj 802 model_ncoeffbound = ii 803 exit 804 end if 805 end if 806 end do 807 end if 808 809 ABI_DEALLOCATE(list_tmp) 810 ABI_DEALLOCATE(list_bound) 811 ABI_DEALLOCATE(isPositive) 812 813 ! Exit if the model is bounded 814 if(effective_potential%anharmonics_terms%bounded) then 815 ! Final transfert 816 write(message, '(3a)' ) ch10,' => The model is now bounded' 817 call wrtout(ab_out,message,'COLL') 818 call wrtout(std_out,message,'COLL') 819 do kk=ncoeff+1,ncoeff+model_ncoeffbound 820 icoeff_bound = listcoeff_bound(model_bound,kk)-ncoeff ! need to remove ncoeff value 821 call polynomial_coeff_free(coeffs_tmp(kk)) 822 call polynomial_coeff_init(coeff_values(model_bound,kk),& 823 & coeffs_bound(icoeff_bound)%nterm,& 824 & coeffs_tmp(kk),& 825 & coeffs_bound(icoeff_bound)%terms,& 826 & coeffs_bound(icoeff_bound)%name,& 827 & check=.false.) 828 end do 829 ABI_DEALLOCATE(coeff_values) 830 ABI_DEALLOCATE(listcoeff_bound) 831 exit 832 end if 833 ABI_DEALLOCATE(coeff_values) 834 835 end do 836 837 do ii=1,ncoeff_bound 838 call polynomial_coeff_free(coeffs_bound(ii)) 839 end do 840 if(allocated(coeffs_bound)) ABI_DEALLOCATE(coeffs_bound) 841 842 end if 843 844 if(.not.effective_potential%anharmonics_terms%bounded)then 845 write(message, '(3a)' ) ch10,' => The model cannot be bounded' 846 call wrtout(ab_out,message,'COLL') 847 call wrtout(std_out,message,'COLL') 848 model_ncoeffbound = 0 849 model_bound = 0 850 end if 851 852 ! Fit the final model 853 call effective_potential_setCoeffs(coeffs_tmp(1:ncoeff+model_ncoeffbound),effective_potential,& 854 & ncoeff+model_ncoeffbound) 855 856 call fit_polynomial_coeff_fit(effective_potential,(/0/),(/0/),hist,0,(/0,0/),0,0,& 857 & -1,1,comm,verbose=.false.) 858 859 write(message, '(3a)') ch10,' Fitted coefficients at the end of the fit bound process: ' 860 call wrtout(ab_out,message,'COLL') 861 call wrtout(std_out,message,'COLL') 862 863 do ii = 1,ncoeff+model_ncoeffbound 864 write(message, '(a,I0,a,ES19.10,2a)') " ",ii," =>",& 865 & effective_potential%anharmonics_terms%coefficients(ii)%coefficient,& 866 & " ",trim(effective_potential%anharmonics_terms%coefficients(ii)%name) 867 call wrtout(ab_out,message,'COLL') 868 call wrtout(std_out,message,'COLL') 869 end do 870 871 ! Deallocation 872 ABI_DEALLOCATE(listcoeff) 873 do ii=1,ncoeff_max 874 call polynomial_coeff_free(coeffs_tmp(ii)) 875 end do 876 if(allocated(coeffs_tmp)) ABI_DEALLOCATE(coeffs_tmp) 877 878 do ii=1,ncoeff_max 879 call polynomial_coeff_free(coeffs_all(ii)) 880 end do 881 if(allocated(coeffs_all)) ABI_DEALLOCATE(coeffs_all) 882 883 end if 884 885 else if (option == -3) then 886 887 !************************************************************* 888 ! Call the routine for calculation of the energy for specific 889 ! partern of displacement or strain for the effective 890 ! Hamiltonian 891 !************************************************************* 892 ! write(message, '((80a),4a)' ) ('-',ii=1,80), ch10,& 893 !& ' Effective Hamiltonian calculation' 894 ! call wrtout(ab_out,message,'COLL') 895 ! call wrtout(std_out,message,'COLL') 896 897 ! acell = one 898 ! call gstate(args_gs,acell,codvsn,cpui,dtfil,dtset,iexit,initialized,& 899 !& mpi_enreg,npwtot,dtset%occ_orig,pawang,pawrad,pawtab,& 900 !& psps,results_gs,dtset%rprimd_orig,scf_history,vel,vel_cell,wvl,xred) 901 902 end if 903 904 !*************************************************************** 905 ! 5 Deallocation of array 906 !*************************************************************** 907 908 ABI_DEALLOCATE(amass) 909 ABI_DEALLOCATE(fred) 910 ABI_DEALLOCATE(fcart) 911 ABI_DEALLOCATE(indsym) 912 ABI_DEALLOCATE(rhog) 913 ABI_DEALLOCATE(rhor) 914 ABI_DEALLOCATE(vel) 915 ABI_DEALLOCATE(xred) 916 ABI_DEALLOCATE(xred_old) 917 ABI_DEALLOCATE(ab_xfh%xfhist) 918 919 ! if(option == -3)then 920 ! call args_gs_free(args_gs) 921 ! call psps_free(psps) 922 ! do ii = 1,npsp 923 ! call paw_setup_free(paw_setup(ii)) 924 ! end do 925 ! ABI_DEALLOCATE(paw_setup) 926 ! ABI_DEALLOCATE(ipsp2xml) 927 ! ABI_DEALLOCATE(pspheads) 928 ! call pawrad_free(pawrad) 929 ! call pawtab_free(pawtab) 930 ! ABI_DATATYPE_DEALLOCATE(pawrad) 931 ! ABI_DATATYPE_DEALLOCATE(pawtab) 932 ! ABI_DEALLOCATE(npwtot) 933 ! end if 934 call dtset_free(dtset) 935 call destroy_results_gs(results_gs) 936 call scfcv_destroy(scfcv_args) 937 call destroy_mpi_enreg(mpi_enreg) 938 939 end if 940 941 write(message, '(a,(80a),a,a)' ) ch10,& 942 & ('=',ii=1,80),ch10 943 call wrtout(ab_out,message,'COLL') 944 call wrtout(std_out,message,'COLL') 945 946 end subroutine mover_effpot