TABLE OF CONTENTS


ABINIT/multibinit [ Programs ]

[ Top ] [ Programs ]

NAME

 multibinit

FUNCTION

 Main routine MULTIBINIT.

COPYRIGHT

 Copyright (C) 1999-2018 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 .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

INPUTS

  (main routine)

OUTPUT

  (main routine)

PARENTS

CHILDREN

      ab7_invars_set_flags,abi_io_redirect,abihist_bcast,abihist_free
      abimem_init,abinit_doctor,compute_anharmonics
      effective_potential_file_getdimsystem,effective_potential_file_gettype
      effective_potential_file_maphisttoref,effective_potential_file_read
      effective_potential_file_readmdfile,effective_potential_free
      effective_potential_setconfinement,effective_potential_writenetcdf
      effective_potential_writexml,fit_polynomial_coeff_fit
      fit_polynomial_printsystemfiles,flush_unit,herald,init10,instrng
      inupper,invars10,isfile,mover_effpot,multibinit_dtset_free
      outvars_multibinit,timein,wrtout,xmpi_bcast,xmpi_init,xmpi_sum

SOURCE

 38 #if defined HAVE_CONFIG_H
 39 #include "config.h"
 40 #endif
 41 
 42 #include "abi_common.h"
 43 
 44 program multibinit
 45 
 46  use defs_basis
 47  use defs_abitypes
 48  use m_build_info
 49  use m_xmpi
 50  use m_xomp
 51  use m_abicore
 52  use m_errors
 53  use m_argparse
 54  use m_effective_potential
 55  use m_fit_polynomial_coeff
 56  use m_multibinit_dataset
 57  use m_effective_potential_file
 58  use m_spin_model
 59  use m_abihist
 60  use m_ab7_invars
 61 
 62  use m_specialmsg, only : specialmsg_getcount, herald
 63  use m_io_tools,   only : flush_unit, open_file
 64  use m_fstrings,   only : replace, inupper
 65  use m_time,       only : asctime, timein
 66  use m_parser,     only : instrng
 67  use m_dtset,      only : chkvars
 68  use m_dtfil,      only : isfile
 69  use m_mover_effpot, only : mover_effpot
 70  !use m_generate_training_set, only : generate_training_set
 71  use m_compute_anharmonics, only : compute_anharmonics
 72  use m_init10,              only : init10
 73 
 74 !This section has been created automatically by the script Abilint (TD).
 75 !Do not modify the following lines by hand.
 76 #undef ABI_FUNC
 77 #define ABI_FUNC 'multibinit'
 78 !End of the abilint section
 79 
 80  implicit none
 81 
 82 !Arguments -----------------------------------
 83 
 84 !Local variables-------------------------------
 85 ! Set array dimensions
 86  integer,parameter :: master=0 ! FIXME: these should not be reserved unit numbers!
 87  integer :: comm,filetype,ii,ierr,lenstr
 88  integer :: natom,nph1l,nrpt,ntypat,nproc,my_rank
 89  integer :: option
 90  logical :: iam_master
 91  real(dp) :: tcpu,tcpui,twall,twalli
 92  real(dp) :: tsec(2)
 93  character(len=24) :: codename,start_datetime
 94  character(len=strlen) :: string
 95  character(len=fnlen) :: filnam(17),tmpfilename,name
 96  character(len=fnlen) :: filstat
 97  character(len=500) :: message
 98  type(multibinit_dtset_type) :: inp
 99  type(effective_potential_type) :: reference_effective_potential
100  type(abihist) :: hist
101  type(args_t) :: args
102 
103 !TODO hexu: add types for spin here.
104  type(spin_model_t) :: spin_model
105 !TEST_AM
106 ! integer :: natom_sp
107 ! real(dp),allocatable :: dynmat(:,:,:,:,:)
108 !TEST_AM
109 !******************************************************************
110 
111 !Change communicator for I/O (mandatory!)
112  call abi_io_redirect(new_io_comm=xmpi_world)
113 
114 !Initialize MPI
115  call xmpi_init()
116 
117 !MPI variables
118  comm = xmpi_world
119  nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
120  iam_master = (my_rank == master)
121 
122 
123 ! Parse command line arguments.
124  args = args_parser(); if (args%exit /= 0) goto 100
125 
126 !Initialize memory profiling if it is activated !if a full abimem.mocc report is desired,
127 !set the argument of abimem_init to "2" instead of "0"
128 !note that abimem.mocc files can easily be multiple GB in size so don't use this option normally
129 #ifdef HAVE_MEM_PROFILING
130  call abimem_init(0)
131 #endif
132 
133 !Initialisation of the timing
134  call timein(tcpui,twalli)
135 
136  if (iam_master) then
137    codename='MULTIBINIT'//repeat(' ',14)
138    call herald(codename,abinit_version,std_out)
139  end if
140 
141  start_datetime = asctime()
142 
143 !Print the number of cpus in log file
144  write(message,'(a,i5,a)') '-  nproc =',nproc,ch10
145  call wrtout(std_out,message,'COLL')
146 
147 !Initialise the code : write heading, and read names of files.
148  call init10(filnam,comm)
149 
150 ! Call the parser from the parser module.
151  filstat = trim("_STATUS")
152  call ab7_invars_set_flags(.true., .true., status_file = filstat, timab_tsec = tsec)
153 
154 !******************************************************************
155 
156  call timein(tcpu,twall)
157 
158  write(message, '(a,f11.3,a,f11.3,a)' )'-begin at tcpu',tcpu-tcpui,'  and twall',twall-twalli,' sec'
159  call wrtout(std_out,message,'COLL')
160 
161 ! Open output files and ab_out (might change its name if needed)
162 ! MJV 1/2010 : now output file is open, but filnam(2) continues unmodified
163 ! so the other output files are overwritten instead of accumulating.
164  if (iam_master) then
165    tmpfilename = filnam(2)
166    call isfile(tmpfilename,'new')
167    if (open_file(tmpfilename,message,unit=ab_out,form="formatted",status="new",&
168 &   action="write") /= 0) then
169      MSG_ERROR(message)
170    end if
171 !  Call open_file(unit=ab_out,file=tmpfilename,form='formatted',status='new')
172    rewind (unit=ab_out)
173    call herald(codename,abinit_version,ab_out)
174 !  Print the number of cpus in output
175    write(message,'(a,i5,a)') '-  nproc =',nproc
176    call wrtout(ab_out,message,'COLL')
177  else
178    ab_out = dev_null
179  end if
180 
181  write(message, '(a,(80a),a)' ) ch10,&
182 & ('=',ii=1,80),ch10
183  call wrtout(ab_out,message,'COLL')
184  call wrtout(std_out,message,'COLL')
185 
186 
187 !To automate a maximum calculation, multibinit reads the number of atoms
188 !in the file (ddb or xml). If DDB file is present in input, the ifc calculation
189 !will be initilaze array to the maximum of atoms (natifc=natom,atifc=1,natom...) in invars10
190  write(message, '(6a)' )' Read the information in the reference structure in ',ch10,&
191 & '-',trim(filnam(3)),ch10,' to initialize the multibinit input'
192  call wrtout(ab_out,message,'COLL')
193  call wrtout(std_out,message,'COLL')
194 
195  call effective_potential_file_getDimSystem(filnam(3),natom,ntypat,nph1l,nrpt)
196 
197 !Read the input file, and store the information in a long string of characters
198 !strlen from defs_basis module
199  option=1
200  if (iam_master) then
201    call instrng (filnam(1),lenstr,option,strlen,string)
202    !To make case-insensitive, map characters to upper case:
203    call inupper(string(1:lenstr))
204 
205    !Check whether the string only contains valid keywords
206    call chkvars(string)
207 
208  end if
209 
210  call xmpi_bcast(string,master, comm, ierr)
211  call xmpi_bcast(lenstr,master, comm, ierr)
212 
213 !Read the input file
214  call invars10(inp,lenstr,natom,string)
215 
216  if (iam_master) then
217 !  Echo the inputs to console and main output file
218    call outvars_multibinit(inp,std_out)
219    call outvars_multibinit(inp,ab_out)
220  end if
221 
222 ! Read and treat the reference structure
223 !****************************************************************************************
224   if (inp%spin_dynamics>0) then
225    if (iam_master) then
226 
227      write(message,'(a,(80a),3a)') ch10,('=',ii=1,80),ch10,ch10,&
228 &     'reading spin terms.'
229      call spin_model_t_initialize(spin_model, filnam(3), inp )
230    end if
231  else
232    !  Read the model (from DDB or XML)  
233    call effective_potential_file_read(filnam(3),reference_effective_potential,inp,comm)
234 
235  !Read the coefficient from fit
236    if(filnam(4)/=''.and.filnam(4)/='no')then
237      call effective_potential_file_getType(filnam(4),filetype)
238    ! TODO hexu: filetype==(33?)
239      if(filetype==3.or.filetype==23) then
240        call effective_potential_file_read(filnam(4),reference_effective_potential,inp,comm)
241      else
242        write(message,'(a,(80a),3a)') ch10,('=',ii=1,80),ch10,ch10,&
243 &         ' There is no specific file for the coefficients from polynomial fitting'
244        call wrtout(ab_out,message,'COLL')
245        call wrtout(std_out,message,'COLL')
246      end if
247    else
248      if(inp%ncoeff/=0) then
249        write(message, '(5a)' )&
250 &         'ncoeff is specified in the input but,',ch10,&
251 &         'there is no file for the coefficients ',ch10,&
252 &         'Action: add coefficients.xml file'
253        MSG_ERROR(message)
254 
255      else
256        write(message,'(a,(80a),3a)') ch10,('=',ii=1,80),ch10,ch10,&
257 &         ' There is no file for the coefficients from polynomial fitting'
258        call wrtout(ab_out,message,'COLL')
259        call wrtout(std_out,message,'COLL')
260      end if
261    end if
262  end if
263  
264 !****************************************************************************************
265 
266 ! Compute the third order derivative with finite differences
267 !****************************************************************************************
268    if (inp%strcpling > 0) then
269      call compute_anharmonics(reference_effective_potential,filnam,inp,comm)
270    end if
271 !****************************************************************************************
272 
273 ! If needed, fit the anharmonic part and compute the confinement potential
274 !****************************************************************************************
275    if (inp%fit_coeff/=0.or.inp%confinement==2.or.inp%bound_model/=0) then
276 
277      if(iam_master) then
278 !    Read the MD file
279        write(message,'(a,(80a),7a)')ch10,('=',ii=1,80),ch10,ch10,&
280 &       '-Reading the file the HIST file :',ch10,&
281 &       '-',trim(filnam(5)),ch10
282 
283        call wrtout(std_out,message,'COLL')
284        call wrtout(ab_out,message,'COLL')
285        if(filnam(5)/=''.and.filnam(5)/='no')then
286          call effective_potential_file_readMDfile(filnam(5),hist,option=inp%ts_option)
287          if (hist%mxhist == 0)then
288            write(message, '(5a)' )&
289 &           'The MD ',trim(filnam(5)),' file is not correct ',ch10,&
290 &           'Action: add MD file'
291            MSG_ERROR(message)
292          end if
293        else
294          if (inp%fit_coeff/=0) then
295            write(message, '(3a)' )&
296 &           'There is no MD file to fit the coefficients ',ch10,&
297 &           'Action: add MD file'
298            MSG_ERROR(message)
299          else
300            if (inp%bound_model/=0) then
301              write(message, '(3a)' )&
302 &             'There is no MD file to bound the model ',ch10,&
303 &             'Action: add MD file'
304              MSG_ERROR(message)
305            else if(inp%confinement==2) then
306              write(message, '(3a)' )&
307 &             'There is no MD file to compute the confinement',ch10,&
308 &             'Action: add MD file'
309              MSG_ERROR(message)
310            end if
311          end if
312        end if
313      end if
314 !  MPI BROADCAST the history of the MD
315      call abihist_bcast(hist,master,comm)
316 !  Map the hist in order to be consistent with the supercell into reference_effective_potential
317      call effective_potential_file_mapHistToRef(reference_effective_potential,hist,comm)
318    end if
319 
320 !TEST_AM
321 ! call effective_potential_checkDEV(reference_effective_potential,hist,size(hist%xred,2),hist%mxhist)
322 ! stop
323 !TEST_AM
324 
325 !Generate the confinement polynome (not working yet)
326    if(inp%confinement/=0)then
327      option=inp%confinement
328      select case(option)
329      case(1)
330        call effective_potential_setConfinement(inp%conf_cutoff_disp,inp%conf_cutoff_strain,&
331 &       reference_effective_potential,inp%conf_power_fact_disp,&
332 &       inp%conf_power_fact_strain,inp%conf_power_disp,&
333 &       inp%conf_power_disp,inp%conf_power_strain,&
334 &       need_confinement=.TRUE.)
335 
336        write(message,'(a,(80a),3a)') ch10,('=',ii=1,80),ch10,ch10,&
337 &       ' The confinement potential is active.'
338        call wrtout(ab_out,message,'COLL')
339        call wrtout(std_out,message,'COLL')
340 
341      case(2)
342        write(message,'(a,(80a),3a)') ch10,('=',ii=1,80),ch10,ch10,&
343 &       ' The confinement potential is computed from the MD file and actived.'
344        call wrtout(ab_out,message,'COLL')
345        call wrtout(std_out,message,'COLL')
346 
347      end select
348    end if
349 
350 
351 !Fit the coeff
352    if (inp%fit_coeff/=0)then
353      option=inp%fit_coeff
354      if(hist%mxhist >0)then
355        if (option==-1)then
356 !      option == -1
357 !      Print the file in the specific format for the script of carlos
358 !      Born_Charges
359 !      Dielectric_Tensor
360 !      harmonic.xml
361 !      Reference_structure
362 !      Strain_Tensor
363 !      symmetry_operations (only cubic)
364          if (iam_master) then
365            call fit_polynomial_printSystemFiles(reference_effective_potential,hist)
366          end if
367        else if (option==1.or.option==2)then
368 !      option = 1
369          call fit_polynomial_coeff_fit(reference_effective_potential,&
370 &         inp%fit_bancoeff,inp%fit_fixcoeff,hist,inp%fit_generateCoeff,&
371 &         inp%fit_rangePower,inp%fit_nbancoeff,inp%fit_ncoeff,&
372 &         inp%fit_nfixcoeff,option,comm,cutoff_in=inp%fit_cutoff,&
373 &         initialize_data=inp%fit_initializeData==1,&
374 &         fit_tolMSDF=inp%fit_tolMSDF,fit_tolMSDS=inp%fit_tolMSDS,fit_tolMSDE=inp%fit_tolMSDE,&
375 &         fit_tolMSDFS=inp%fit_tolMSDFS,&
376 &         verbose=.true.,positive=.false.,&
377 &         anharmstr=inp%fit_anhaStrain==1,&
378 &         spcoupling=inp%fit_SPCoupling==1)
379        end if
380      else
381        write(message, '(3a)' )&
382 &       'There is no step in the MD file ',ch10,&
383 &       'Action: add correct MD file'
384        MSG_ERROR(message)
385      end if
386    end if
387 
388 
389 !try to bound the model with mover_effpot
390 !we need to use the molecular dynamics
391    if(inp%bound_model>0.and.inp%bound_model<=2)then
392      call mover_effpot(inp,filnam,reference_effective_potential,-1*inp%bound_model,comm,hist=hist)
393    end if
394 
395 !****************************************************************************************
396 
397 !****************************************************************************************
398 !TEST_AM
399 !Effective Hamiltonian, compute the energy for given patern
400 ! call mover_effpot(inp,filnam,reference_effective_potential,-2,comm,hist=hist)
401 !TEST_AM
402 
403 !****************************************************************************************
404 
405 !****************************************************************************************
406 !Print the effective potential system + coefficients (only master CPU)
407 ! TODO hexu: add print spin model.
408    if(iam_master) then
409      if (inp%prt_model >= 1) then
410        write(message, '(a,(80a),a)' ) ch10,&
411 &       ('=',ii=1,80)
412        call wrtout(ab_out,message,'COLL')
413        call wrtout(std_out,message,'COLL')
414        name = replace(trim(filnam(2)),".out","")
415        call effective_potential_writeXML(reference_effective_potential,inp%prt_model,filename=name,&
416 &       prt_dipdip=inp%dipdip_prt==1)
417      else if (inp%prt_model == -2)then
418 !    NetCDF case, in progress
419        name = trim(filnam(2))//"_sys.nc"
420        call effective_potential_writeNETCDF(reference_effective_potential,1,filename=name)
421      end if
422    end if
423 !****************************************************************************************
424 
425 !TEST_AM SECTION
426 ! Print the Phonon dos/spectrum
427 ! if(inp%prt_phfrq > 0) then
428 !     call effective_potential_printPDOS(reference_effective_potential,filnam(2),&
429 !&           inp%ncell,inp%nph1l,inp%prt_phfrq,inp%qph1l)
430 !   end if
431 
432 !Intialisation of the effective potential type
433 !   call effective_potential_file_read(filnam(4),reference_effective_potential,inp,comm)
434 !   name = "test.xml"
435 !   call effective_potential_writeXML(reference_effective_potential,1,filename=name)
436 ! just for TEST
437 !   if(inp%prt_phfrq > 0) then
438 !      natom_sp = reference_effective_potential%supercell%natom_supercell
439 !      ABI_ALLOCATE(dynmat,(2,3,natom_sp,3,natom_sp))
440 !      call effective_potential_effpot2dynmat(dynmat,inp%delta_df,reference_effective_potential,&
441 ! &                                           reference_effective_potential%supercell%natom_supercell,&
442 ! &                                           int(reference_effective_potential%supercell%qphon),3)
443 
444 !      ABI_DEALLOCATE(dynmat)
445 !    end if
446 ! end if
447 !TEST_AM SECTION
448 
449 
450 ! Compute the monte carlo, molecular dynamics of compute specific energy
451 !****************************************************************************************
452    if(inp%dynamics>=1) then
453      call mover_effpot(inp,filnam,reference_effective_potential,inp%dynamics,comm)
454    end if
455 
456 !****************************************************************************************
457 
458 
459 ! Run spin dynamics
460 !****************************************************************************************
461    if(inp%spin_dynamics>0) then
462   ! TODO hexu: no mpi yet.
463      if(iam_master) then
464        call spin_model_t_run(spin_model)
465      end if
466    end if
467 !****************************************************************************************
468 
469 
470 !Free the effective_potential and dataset
471 !****************************************************************************************
472    call effective_potential_free(reference_effective_potential)
473    call multibinit_dtset_free(inp)
474    call abihist_free(hist)
475    call spin_model_t_finalize(spin_model)
476 !****************************************************************************************
477 
478    write(message,'(a,a,a,(80a))') ch10,('=',ii=1,80),ch10
479    call wrtout(ab_out,message,'COLL')
480    call wrtout(std_out,message,'COLL')
481 
482    call timein(tcpu,twall)
483    tsec(1)=tcpu-tcpui
484    tsec(2)=twall-twalli
485 
486    write(message, '(a,i4,a,f13.1,a,f13.1)' )' Proc.',my_rank,' individual time (sec): cpu=',&
487 &   tsec(1),'  wall=',tsec(2)
488    call wrtout(std_out,message,"COLL")
489 
490    if (iam_master) then
491      write(ab_out, '(a,a,a,i4,a,f13.1,a,f13.1)' )'-',ch10,&
492 &     '- Proc.',my_rank,' individual time (sec): cpu=',tsec(1),'  wall=',tsec(2)
493    end if
494 
495    call xmpi_sum(tsec,comm,ierr)
496 
497    write(message, '(a,(80a),a,a,a,f11.3,a,f11.3,a,a,a,a)' ) ch10,&
498 &   ('=',ii=1,80),ch10,ch10,&
499 &   '+Total cpu time',tsec(1),&
500 &   '  and wall time',tsec(2),' sec',ch10,ch10,&
501 &   ' multibinit : the run completed succesfully.'
502    call wrtout(std_out,message,'COLL')
503    call wrtout(ab_out,message,'COLL')
504 
505    if (iam_master) then
506    ! Write YAML document with the final summary.
507    ! we use this doc to test whether the calculation is completed.
508      write(std_out,"(a)")"--- !FinalSummary"
509      write(std_out,"(a)")"program: multibinit"
510      write(std_out,"(2a)")"version: ",trim(abinit_version)
511      write(std_out,"(2a)")"start_datetime: ",start_datetime
512      write(std_out,"(2a)")"end_datetime: ",asctime()
513      write(std_out,"(a,f13.1)")"overall_cpu_time: ",tsec(1)
514      write(std_out,"(a,f13.1)")"overall_wall_time: ",tsec(2)
515      write(std_out,"(a,i0)")"mpi_procs: ",xmpi_comm_size(xmpi_world)
516      write(std_out,"(a,i0)")"omp_threads: ",xomp_get_num_threads(open_parallel=.True.)
517    !write(std_out,"(a,i0)")"num_warnings: ",nwarning
518    !write(std_out,"(a,i0)")"num_comments: ",ncomment
519      write(std_out,"(a)")"..."
520      call flush_unit(std_out)
521    end if
522 
523 !Write information on file about the memory before ending mpi module, if memory profiling is enabled
524    call abinit_doctor("__multibinit")
525 
526    call flush_unit(ab_out)
527    call flush_unit(std_out)
528 
529    if (iam_master) close(ab_out)
530 
531    100 call xmpi_end()
532 
533    end program multibinit