Line data Source code
1 :
2 : !!****m* ABINIT/m_multibinit_main
3 : !! NAME
4 : !! m_multibinit_main
5 : !!
6 : !! FUNCTION
7 : !! Main routine MULTIBINIT.
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 1999-2026 ABINIT group (AM, hexu)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
15 : !!
16 : !! INPUTS
17 : !! (main routine)
18 : !!
19 : !! OUTPUT
20 : !! (main routine)
21 : !! NOTES
22 : !! Should be
23 : !! 1 moved to somewhere else
24 : !! 2 be replaced with the new implementation multibinit_main2.
25 : !!
26 : !! SOURCE
27 :
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
72 : !!***
73 : contains
74 : !!****f* m_multbinit_main/multibinit_main
75 : !!
76 : !! NAME
77 : !! multibinit_main
78 : !!
79 : !! FUNCTION
80 : !! The main function of multibinit
81 : !!
82 : !! INPUTS
83 : !! filnam: The filenames from the files file. 17 files in total.
84 : !!
85 : !! OUTPUT
86 : !!
87 : !! SOURCE
88 67 : 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 67 : type(multibinit_dtset_type), target :: inp
93 7571 : type(effective_potential_type) :: reference_effective_potential, read_effective_potential
94 67 : 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 67 : master = 0
143 67 : comm = xmpi_world
144 67 : nproc = xmpi_comm_size(comm)
145 67 : my_rank = xmpi_comm_rank(comm)
146 67 : 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 67 : option=1
151 67 : if (iam_master) then
152 43 : call instrng (filnam(1),lenstr,option,strlen,string,raw_string)
153 : !To make case-insensitive, map characters to upper case:
154 43 : call inupper(string(1:lenstr))
155 : !Check whether the string only contains valid keywords
156 43 : call chkvars(string)
157 : end if
158 67 : call xmpi_bcast(string,master, comm, ierr)
159 67 : call xmpi_bcast(raw_string,master, comm, ierr)
160 67 : 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 67 : if(trim(filnam(3)) /='') then
170 67 : sys_fname=filnam(3)
171 : else
172 0 : 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 67 : if (iam_master) then
179 43 : write(message, '(6a)' )' Read the information in the reference structure in ',ch10,&
180 86 : & '-',trim(sys_fname),ch10,' to initialize the multibinit input'
181 43 : call wrtout(ab_out,message,'COLL')
182 43 : call wrtout(std_out,message,'COLL')
183 : end if
184 :
185 67 : 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 67 : call invars10(inp,lenstr,natom,string)
191 67 : call postfix_fnames(input_path, filnam, inp)
192 :
193 67 : need_new_multibinit= inp%spin_dynamics > 0 .or. inp%lwf_dynamics > 0 .or. inp%dynamics >= 100
194 :
195 67 : if (iam_master) then
196 43 : if(need_new_multibinit) then
197 0 : 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 43 : call outvars_multibinit(inp,std_out)
201 43 : call outvars_multibinit(inp,ab_out)
202 : end if
203 :
204 67 : if(dry_run/=0) then
205 0 : call wrtout([std_out, ab_out], "Multibinit in dry_run mode. Exiting after input parser")
206 0 : call xmpi_end()
207 : !goto 100
208 : endif
209 :
210 : ! Read the model (from DDB or XML)
211 67 : call effective_potential_file_read(filnam(3),reference_effective_potential,inp,comm)
212 : ! TODO: coeff_file_rw is currently set to 0 as this functionality needs more testing.
213 67 : if (inp%coeff_file_rw==2) then
214 0 : read_effective_potential=reference_effective_potential
215 : !call effective_potential_file_read(filnam(3),read_effective_potential,inp,comm)
216 0 : if(filnam(6)/=''.and.filnam(6)/='no')then
217 0 : call effective_potential_file_getType(filnam(6),filetype)
218 0 : if(filetype==3.or.filetype==23) then
219 0 : call effective_potential_file_read(filnam(6),read_effective_potential,inp,comm)
220 : end if
221 : end if
222 : end if
223 :
224 67 : if(filnam(4)/=''.and.filnam(4)/='no') then
225 24 : call effective_potential_file_getType(filnam(4),filetype)
226 24 : if(filetype==3.or.filetype==23) then
227 24 : call effective_potential_file_read(filnam(4),reference_effective_potential,inp,comm)
228 : else
229 0 : write(message,'(a,(80a),3a)') ch10,('=',ii=1,80),ch10,ch10,&
230 0 : & ' There is no specific file for the coefficients from polynomial fitting'
231 0 : call wrtout(ab_out,message,'COLL')
232 0 : call wrtout(std_out,message,'COLL')
233 : end if
234 : else
235 43 : if(inp%ncoeff/=0) then
236 : write(message, '(5a)' )&
237 0 : & 'ncoeff is specified in the input but,',ch10,&
238 0 : & 'there is no file for the coefficients ',ch10,&
239 0 : & 'Action: add coefficients.xml file'
240 0 : ABI_ERROR(message)
241 :
242 : else
243 3526 : write(message,'(a,(80a),3a)') ch10,('=',ii=1,80),ch10,ch10,&
244 3569 : & ' There is no file for the coefficients from polynomial fitting'
245 43 : call wrtout(ab_out,message,'COLL')
246 43 : call wrtout(std_out,message,'COLL')
247 : end if
248 : end if
249 :
250 : !****************************************************************************************
251 : !SCALE UP Initialize the electronic model (If scale-up is available)
252 : !****************************************************************************************
253 67 : elec_eval = .FALSE.
254 :
255 : #if defined DEV_MS_SCALEUP
256 : if(inp%scup_dtset%scup_elec_model)then
257 : write(message,'(a,(80a),4a)') ch10,('=',ii=1,80),ch10,ch10,&
258 : ' Initializing Electronic Model with SCALE-UP',ch10
259 : call wrtout(ab_out,message,'COLL')
260 : call wrtout(std_out,message,'COLL')
261 :
262 : !Set Variables
263 : elec_eval = .TRUE.
264 : ksamp = inp%scup_dtset%scup_ksamp
265 : tcharge = inp%scup_dtset%scup_tcharge
266 : if(inp%scup_dtset%scup_ismagnetic)ismagnetic=.TRUE.
267 : if(inp%scup_dtset%scup_istddft)istddft=.TRUE.
268 : if(inp%scup_dtset%scup_initorbocc)initorbocc=.TRUE.
269 :
270 : ! Call to Scale-Up
271 : err_init_elec = global_init_model(filnam(3),inp%ncell,needlattice,needelectrons,didi,&
272 : & harm_der,tcharge,ksamp,ismagnetic,istddft,initorbocc)
273 :
274 : !Set Print variables
275 : if(inp%scup_dtset%scup_printgeom)printgeom=.TRUE.
276 : if(inp%scup_dtset%scup_printeigv)printeigv=.TRUE.
277 : if(inp%scup_dtset%scup_printeltic)printeltic=.TRUE.
278 : if(inp%scup_dtset%scup_printorbocc)printorbocc=.TRUE.
279 :
280 : !Set Print Parameters within scaleup
281 : call global_set_print_parameters(printgeom,printeigv,printeltic,printorbocc,&
282 : & inp%scup_dtset%scup_printbands)
283 :
284 : !Set SCF controling variables (values contain defaults, if not specified in the input)
285 : call global_set_scf_parameters(inp%scup_dtset%scup_scfmixing,inp%scup_dtset%scup_scfthresh,&
286 : & inp%scup_dtset%scup_smearing,inp%scup_dtset%scup_maxscfstep,&
287 : & inp%scup_dtset%scup_startpulay,inp%scup_dtset%scup_freezden)
288 :
289 :
290 : !Create kpath if printbands=true and pass it to SCALE UP
291 : if(inp%scup_dtset%scup_printbands)then
292 :
293 : call scup_kpath_new(inp%scup_dtset%scup_speck,&
294 : & reference_effective_potential%supercell%rprimd,&
295 : & inp%scup_dtset%scup_ndivsm,inp%scup_dtset%scup_kpath)
296 : call scup_kpath_print(inp%scup_dtset%scup_kpath)
297 :
298 : call global_set_print_bands(inp%scup_dtset%scup_printbands,&
299 : & inp%scup_dtset%scup_nspeck,inp%scup_dtset%scup_kpath%ndivs,&
300 : & inp%scup_dtset%scup_speck)
301 : endif
302 : endif
303 : #endif
304 :
305 : !****************************************************************************************
306 : ! Compute the third order derivative with finite differences
307 : !****************************************************************************************
308 67 : if (inp%strcpling > 0) then
309 0 : call compute_anharmonics(reference_effective_potential,filnam,inp,comm)
310 : end if
311 : !****************************************************************************************
312 :
313 : ! If needed, fit the anharmonic part and compute the confinement potential
314 : !****************************************************************************************
315 67 : if (inp%fit_coeff/=0.or.inp%confinement==2.or.inp%bound_model/=0 .or. inp%opt_effpot/=0) then
316 :
317 20 : if(iam_master) then
318 : ! Read the MD file
319 984 : write(message,'(a,(80a),7a)')ch10,('=',ii=1,80),ch10,ch10,&
320 12 : & '-Reading the training-set file :',ch10,&
321 996 : & '-',trim(filnam(5)),ch10
322 :
323 12 : call wrtout(std_out,message,'COLL')
324 12 : call wrtout(ab_out,message,'COLL')
325 12 : if(filnam(5)/=''.and.filnam(5)/='no')then
326 12 : call effective_potential_file_readMDfile(filnam(5),hist,option=inp%ts_option)
327 :
328 12 : if (hist%mxhist == 0)then
329 : write(message, '(5a)' )&
330 0 : & 'The trainig-set ',trim(filnam(5)),' file is not correct ',ch10,&
331 0 : & 'Action: add training-set file'
332 0 : ABI_ERROR(message)
333 : end if
334 : else
335 0 : if (inp%fit_coeff/=0) then
336 : write(message, '(3a)' )&
337 0 : & 'There is no training-set file to fit the lattice model ',ch10,&
338 0 : & 'Action: add trainings-set-file'
339 0 : ABI_ERROR(message)
340 0 : else if (inp%bound_model/=0) then
341 : write(message, '(3a)' )&
342 0 : & 'There is no training-set file to bound the model ',ch10,&
343 0 : & 'Action: add training-set file '
344 0 : ABI_ERROR(message)
345 0 : else if(inp%confinement==2) then
346 : write(message, '(3a)' )&
347 0 : & 'There is no training-set file to compute the confinement',ch10,&
348 0 : & 'Action: add training-set file '
349 0 : ABI_ERROR(message)
350 0 : else if(inp%opt_effpot==2) then
351 : write(message, '(3a)' )&
352 0 : & 'There is no training-set file to optimize the latice model',ch10,&
353 0 : & 'Action: add training-set file '
354 0 : ABI_ERROR(message)
355 : end if
356 : end if
357 : end if
358 : ! MPI BROADCAST the history of the MD
359 20 : call abihist_bcast(hist,master,comm)
360 : ! Map the hist in order to be consistent with the supercell into reference_effective_potential
361 20 : call effective_potential_file_mapHistToRef(reference_effective_potential,hist,comm)
362 :
363 : end if
364 :
365 : !TEST_AM
366 : ! call effective_potential_checkDEV(reference_effective_potential,hist,size(hist%xred,2),hist%mxhist)
367 : ! stop
368 : !TEST_AM
369 :
370 : !Generate the confinement polynome (not working yet)
371 67 : if(inp%confinement/=0)then
372 0 : option=inp%confinement
373 0 : select case(option)
374 : case(1)
375 : call effective_potential_setConfinement(inp%conf_cutoff_disp,inp%conf_cutoff_strain,&
376 : & reference_effective_potential,inp%conf_power_fact_disp,&
377 : & inp%conf_power_fact_strain,inp%conf_power_disp,&
378 : & inp%conf_power_disp,inp%conf_power_strain,&
379 0 : & need_confinement=.TRUE.)
380 :
381 0 : write(message,'(a,(80a),3a)') ch10,('=',ii=1,80),ch10,ch10,&
382 0 : & ' The confinement potential is active.'
383 0 : call wrtout(ab_out,message,'COLL')
384 0 : call wrtout(std_out,message,'COLL')
385 :
386 : case(2)
387 0 : write(message,'(a,(80a),3a)') ch10,('=',ii=1,80),ch10,ch10,&
388 0 : & ' The confinement potential is computed from the MD file and actived.'
389 0 : call wrtout(ab_out,message,'COLL')
390 0 : call wrtout(std_out,message,'COLL')
391 :
392 : end select
393 : end if
394 :
395 :
396 : !Fit the coeff
397 67 : if (inp%fit_coeff/=0)then
398 18 : option=inp%fit_coeff
399 18 : if(hist%mxhist >0)then
400 18 : if (option==-1)then
401 : ! option == -1
402 : ! Print the file in the specific format for the script of carlos
403 : ! Born_Charges
404 : ! Dielectric_Tensor
405 : ! harmonic.xml
406 : ! Reference_structure
407 : ! Strain_Tensor
408 : ! symmetry_operations (only cubic)
409 0 : if (iam_master) then
410 0 : call fit_polynomial_printSystemFiles(reference_effective_potential,hist)
411 : end if
412 18 : else if (option==1.or.option==2)then
413 : ! option = 1
414 18 : if(inp%fit_iatom/=0)then
415 : call fit_polynomial_coeff_fit(reference_effective_potential,&
416 : & inp%fit_bancoeff,inp%fit_fixcoeff,hist,inp%fit_generateCoeff,&
417 : & inp%fit_rangePower,inp%fit_nbancoeff,inp%fit_ncoeff,&
418 : & inp%fit_nfixcoeff,inp%fit_nimposecoeff,inp%fit_imposecoeff,&
419 : & option,comm,cutoff_in=inp%fit_cutoff,&
420 : & max_power_strain=inp%fit_SPC_maxS,initialize_data=inp%fit_initializeData==1,&
421 : & fit_tolMSDF=inp%fit_tolMSDF,fit_tolMSDS=inp%fit_tolMSDS,fit_tolMSDE=inp%fit_tolMSDE,&
422 : & fit_tolMSDFS=inp%fit_tolMSDFS,fit_tolGF=inp%fit_tolGF,&
423 : & verbose=.true.,positive=.false.,&
424 : & anharmstr=inp%fit_anhaStrain==1,&
425 : & spcoupling=inp%fit_SPCoupling==1,prt_anh=inp%analyze_anh_pot,&
426 : & fit_iatom=inp%fit_iatom,prt_files=.TRUE.,fit_on=inp%fit_on,sel_on=inp%sel_on,&
427 : & fit_factors=inp%fit_factors,prt_GF_csv=inp%prt_GF_csv,dispterms=inp%fit_dispterms==1,&
428 : & coeff_file_rw=inp%coeff_file_rw,read_effective_potential=read_effective_potential, &
429 : & max_nbody=inp%fit_max_nbody, &
430 8 : & drop_rate=inp%fit_drop_rate, ncoeff_per_cycle=inp%fit_ncoeff_per_cycle, fit_weight_T=inp%fit_weight_T)
431 : else
432 10 : if (inp%fit_ncoeff_per_iatom/=0)then
433 1 : if (mod(inp%fit_ncoeff,inp%fit_ncoeff_per_iatom) /= 0)then
434 0 : write(message,'(2a,I3,2a,I3,3a)') ch10,&
435 0 : & 'fit_ncoeff_per_iatom = ', inp%fit_ncoeff_per_iatom,ch10,&
436 0 : & 'is not a divider of fit_ncoeff = ', inp%fit_ncoeff,ch10,&
437 0 : & 'Action: Change fit_ncoeff and/or fit_ncoeff_per_iatom',ch10
438 0 : ABI_ERROR(message)
439 : endif
440 1 : niter = inp%fit_ncoeff/inp%fit_ncoeff_per_iatom
441 1 : if (mod(niter,reference_effective_potential%crystal%nirredat) /= 0)then
442 0 : write(message,'(2a,I3,2a,I3,2a,I3,3a)') ch10,&
443 0 : & 'fit_ncoeff_per_iatom = ', inp%fit_ncoeff_per_iatom,ch10,&
444 0 : & 'times the number of irreducible atoms = ',reference_effective_potential%crystal%nirredat,ch10,&
445 0 : & 'is not a divider of fit_ncoeff = ', inp%fit_ncoeff,ch10,&
446 0 : & 'Action: Change fit_ncoeff and/or fit_ncoeff_per_iatom',ch10
447 0 : ABI_ERROR(message)
448 : endif
449 1 : niter = niter/reference_effective_potential%crystal%nirredat
450 : else if (inp%fit_ncoeff_per_iatom == 0)then
451 9 : if (mod(inp%fit_ncoeff,reference_effective_potential%crystal%nirredat) /= 0)then
452 0 : write(message,'(2a,I3,2a,I3,3a)') ch10,&
453 0 : & 'The number of irreducible atoms = ',reference_effective_potential%crystal%nirredat,ch10,&
454 0 : & 'is not a divider of fit_ncoeff = ', inp%fit_ncoeff,ch10,&
455 0 : & 'Action: Change fit_ncoeff',ch10
456 0 : ABI_ERROR(message)
457 : endif
458 9 : inp%fit_ncoeff_per_iatom = inp%fit_ncoeff/reference_effective_potential%crystal%nirredat
459 9 : niter = 1
460 : endif
461 820 : write(message,'(a,(80a),7a,I3,3a,I3,3a,I3,3a,I3,2a)') ch10,('=',ii=1,80),ch10,ch10,&
462 10 : & ' Starting Fit Iterations ',ch10,&
463 10 : & ' ----------------------- ',ch10,&
464 10 : & ' Select in total fit_ncoeff = ', inp%fit_ncoeff,' coefficients',ch10,&
465 10 : & ' In ', niter,' iterations',ch10,&
466 10 : & ' Over ', reference_effective_potential%crystal%nirredat, ' irreducible atoms',ch10,&
467 830 : & ' Selecting ', inp%fit_ncoeff_per_iatom, ' coefficients per atom in each iteration',ch10
468 10 : call wrtout(std_out,message,'COLL')
469 10 : call wrtout(ab_out,message,'COLL')
470 10 : need_prt_files=.FALSE.
471 20 : do iiter=1,niter
472 820 : write(message,'(a,(80a),3a,I3,a,I3,2a)') ch10,('-',ii=1,80),ch10,ch10,&
473 830 : & ' Start Iteration (',iiter,'/',niter,')',ch10
474 10 : call wrtout(std_out,message,'COLL')
475 10 : call wrtout(ab_out,message,'COLL')
476 50 : do ii=1,reference_effective_potential%crystal%nirredat
477 30 : if(ii == reference_effective_potential%crystal%nirredat .and. iiter==niter)need_prt_files=.TRUE.
478 30 : if(ii > 1 .or. iiter > 1)inp%fit_nfixcoeff = -1
479 : call fit_polynomial_coeff_fit(reference_effective_potential,&
480 : & inp%fit_bancoeff,inp%fit_fixcoeff,hist,inp%fit_generateCoeff,&
481 : & inp%fit_rangePower,inp%fit_nbancoeff,inp%fit_ncoeff_per_iatom,&
482 : & inp%fit_nfixcoeff,inp%fit_nimposecoeff,inp%fit_imposecoeff,&
483 : & option,comm,cutoff_in=inp%fit_cutoff,&
484 : & max_power_strain=inp%fit_SPC_maxS,initialize_data=inp%fit_initializeData==1,&
485 : & fit_tolMSDF=inp%fit_tolMSDF,fit_tolMSDS=inp%fit_tolMSDS,fit_tolMSDE=inp%fit_tolMSDE,&
486 : & fit_tolMSDFS=inp%fit_tolMSDFS,fit_tolGF=inp%fit_tolGF,&
487 : & verbose=.true.,positive=.false.,&
488 : & anharmstr=inp%fit_anhaStrain==1,&
489 : & spcoupling=inp%fit_SPCoupling==1,prt_anh=inp%analyze_anh_pot,&
490 : & fit_iatom=reference_effective_potential%crystal%irredatindx(ii),&
491 : & prt_files=need_prt_files,fit_on=inp%fit_on,sel_on=inp%sel_on,&
492 : & fit_factors=inp%fit_factors,prt_GF_csv=inp%prt_GF_csv,dispterms=inp%fit_dispterms==1, &
493 : & coeff_file_rw=inp%coeff_file_rw,read_effective_potential=read_effective_potential, &
494 : & max_nbody=inp%fit_max_nbody, &
495 40 : & drop_rate=inp%fit_drop_rate, ncoeff_per_cycle=inp%fit_ncoeff_per_cycle, fit_weight_T=inp%fit_weight_T)
496 : enddo
497 : enddo
498 : endif
499 : end if
500 : else
501 : write(message, '(3a)' )&
502 0 : & 'There is no step in the MD file ',ch10,&
503 0 : & 'Action: add correct MD file'
504 0 : ABI_ERROR(message)
505 : end if
506 : end if
507 :
508 :
509 : !try to bound the model with mover_effpot
510 : !we need to use the molecular dynamics
511 67 : if(inp%bound_model>0.and.inp%bound_model<=2)then
512 0 : call mover_effpot(inp,filnam,reference_effective_potential,-1*inp%bound_model,comm,hist=hist)
513 : !Marcus: New option for bound_model: use optimize routine for generting specific high order terms
514 67 : elseif(inp%bound_model == 3)then
515 82 : write(message,'(a,(80a),4a)')ch10,('=',ii=1,80),ch10,ch10,&
516 83 : & 'Bound Process 3: Generate equivalent high order terms',ch10
517 1 : call wrtout(std_out,message,'COLL')
518 1 : call wrtout(ab_out,message,'COLL')
519 :
520 : call opt_effpotbound(reference_effective_potential,inp%bound_rangePower,hist, inp%bound_EFS,&
521 1 : & inp%bound_factors,inp%bound_penalty,comm, fit_weight_T=inp%fit_weight_T, bound_option=inp%bound_option)
522 :
523 : end if
524 :
525 : !****************************************************************************************
526 : ! OPTIMIZE SECTION, Optimize selected coefficients of effective potential while
527 : ! keeping the others constant
528 : !****************************************************************************************
529 :
530 67 : if(inp%opt_effpot == 1)then
531 82 : write(message,'(a,(80a),4a)')ch10,('=',ii=1,80),ch10,ch10,&
532 83 : & 'Optimizing Effective Potential',ch10
533 :
534 1 : call wrtout(std_out,message,'COLL')
535 1 : call wrtout(ab_out,message,'COLL')
536 :
537 1 : need_analyze_anh_pot = .FALSE.
538 1 : if(inp%analyze_anh_pot == 1) need_analyze_anh_pot = .TRUE.
539 :
540 : call opt_effpot(reference_effective_potential,inp%opt_ncoeff,inp%opt_coeff,hist,inp%opt_on,&
541 1 : & inp%opt_factors,comm,fit_weight_T=inp%fit_weight_T, print_anh=need_analyze_anh_pot)
542 : end if
543 :
544 :
545 :
546 : !****************************************************************************************
547 : ! TEST SECTION test effective potential with regard to test-set
548 : !****************************************************************************************
549 67 : if(inp%test_effpot == 1)then
550 1 : if(iam_master) then
551 : ! Read the test-set .nc file
552 82 : write(message,'(a,(80a),9a)')ch10,('=',ii=1,80),ch10,ch10,&
553 1 : & 'TEST - SET Option',ch10,&
554 1 : & '-Reading the test-set file :',ch10,&
555 83 : & '-',trim(filnam(6)),ch10
556 :
557 1 : call wrtout(std_out,message,'COLL')
558 1 : call wrtout(ab_out,message,'COLL')
559 1 : if(filnam(6)/=''.and.filnam(6)/='no')then
560 1 : call effective_potential_file_readMDfile(filnam(6),hist_tes,option=inp%ts_option)
561 1 : if (hist_tes%mxhist == 0)then
562 : write(message, '(5a)' )&
563 0 : & 'The test-set ',trim(filnam(6)),' file is empty ',ch10,&
564 0 : & 'Action: add non-empty test-set'
565 0 : ABI_ERROR(message)
566 : end if
567 : else
568 : write(message, '(3a)' )&
569 0 : & 'There is no test-set file ',ch10,&
570 0 : & 'Action: add test-set file'
571 0 : ABI_ERROR(message)
572 : end if
573 : end if
574 : ! MPI BROADCAST the history of the MD
575 1 : call abihist_bcast(hist_tes,master,comm)
576 : ! Map the hist in order to be consistent with the supercell into reference_effective_potential
577 1 : call effective_potential_file_mapHistToRef(reference_effective_potential,hist_tes,comm)
578 : ! Initialize if to print anharmonic contribution to energy or not
579 1 : need_analyze_anh_pot = .FALSE.
580 1 : if(inp%analyze_anh_pot == 1) need_analyze_anh_pot = .TRUE.
581 : ! Call to test routine
582 : call fit_polynomial_coeff_testEffPot(reference_effective_potential,hist_tes,master,comm,&
583 : & print_anharmonic=need_analyze_anh_pot,scup_dtset=inp%scup_dtset,&
584 1 : & prt_ph=inp%test_prt_ph)
585 :
586 :
587 :
588 : end if ! End if(inp%test_effpot == 1)then
589 :
590 : !TEST_AM
591 : !Effective Hamiltonian, compute the energy for given patern
592 : ! call mover_effpot(inp,filnam,reference_effective_potential,-2,comm,hist=hist)
593 : !TEST_AM
594 :
595 : !****************************************************************************************
596 :
597 : !****************************************************************************************
598 : !Print the effective potential system + coefficients (only master CPU)
599 67 : if(iam_master) then
600 43 : if (inp%prt_model >= 1) then
601 29 : write(message, '(a,(80a),a)' ) ch10,&
602 2349 : & ('=',ii=1,80)
603 29 : call wrtout(ab_out,message,'COLL')
604 29 : call wrtout(std_out,message,'COLL')
605 : !name = replace(trim(filnam(2)),".out","")
606 : ! Assume new .abo convention
607 29 : name = replace(trim(filnam(2)),".abo","")
608 : call effective_potential_writeXML(reference_effective_potential,inp%prt_model,filename=name,&
609 29 : & prt_dipdip=inp%dipdip_prt==1)
610 14 : else if (inp%prt_model == -2)then
611 : ! NetCDF case, in progress
612 0 : name = trim(filnam(2))//"_sys.nc"
613 0 : call effective_potential_writeNETCDF(reference_effective_potential,1,filename=name)
614 : end if
615 : end if
616 : !****************************************************************************************
617 :
618 : !TEST_AM SECTION
619 : ! Print the Phonon dos/spectrum
620 : ! if(inp%prt_phfrq > 0) then
621 : ! call effective_potential_printPDOS(reference_effective_potential,filnam(2),&
622 : !& inp%ncell,inp%nph1l,inp%prt_phfrq,inp%qph1l)
623 : ! end if
624 :
625 : !Intialisation of the effective potential type
626 : ! call effective_potential_file_read(filnam(4),reference_effective_potential,inp,comm)
627 : ! name = "test.xml"
628 : ! call effective_potential_writeXML(reference_effective_potential,1,filename=name)
629 : ! just for TEST
630 : ! if(inp%prt_phfrq > 0) then
631 : ! natom_sp = reference_effective_potential%supercell%natom_supercell
632 : ! ABI_MALLOC(dynmat,(2,3,natom_sp,3,natom_sp))
633 : ! call effective_potential_effpot2dynmat(dynmat,inp%delta_df,reference_effective_potential,&
634 : ! & reference_effective_potential%supercell%natom_supercell,&
635 : ! & int(reference_effective_potential%supercell%qphon),3)
636 :
637 : ! ABI_FREE(dynmat)
638 : ! end if
639 : ! end if
640 : !TEST_AM SECTION
641 :
642 :
643 : ! Run lattice dynamics (relaxation or molecular dynamics, most of abinits ionmovs are allowed)
644 : !****************************************************************************************
645 67 : if(inp%dynamics>=1) then
646 44 : call mover_effpot(inp,filnam,reference_effective_potential,inp%dynamics,comm)
647 : end if
648 :
649 : !****************************************************************************************
650 :
651 :
652 :
653 : !Free the effective_potential and dataset
654 : !****************************************************************************************
655 67 : call effective_potential_free(reference_effective_potential)
656 67 : call multibinit_dtset_free(inp)
657 67 : call abihist_free(hist)
658 67 : call abihist_free(hist_tes)
659 : !****************************************************************************************
660 :
661 67 : end subroutine multibinit_main
662 : !!***
663 :
664 : end module m_multibinit_driver
|
