Branch data 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
|
