Line data Source code
1 : !!****m* ABINIT/m_mover_effpot
2 : !! NAME
3 : !! m_mover_effpot
4 : !!
5 : !! FUNCTION
6 : !!
7 : !! COPYRIGHT
8 : !! Copyright (C) 2008-2026 ABINIT group (AM)
9 : !! This file is distributed under the terms of the
10 : !! GNU General Public License, see ~abinit/COPYING
11 : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : !!
13 : !! SOURCE
14 :
15 : #if defined HAVE_CONFIG_H
16 : #include "config.h"
17 : #endif
18 :
19 : #include "abi_common.h"
20 :
21 : module m_mover_effpot
22 :
23 : use defs_basis
24 : use m_errors
25 : use m_abicore
26 : use m_dtset
27 : use m_dtfil
28 : use m_abimover
29 : use m_scf_history
30 : use defs_wvltypes
31 : use m_xmpi
32 : use m_phonons
33 : use m_strain
34 : use m_effective_potential_file
35 : use m_supercell
36 : use m_psps
37 : use m_args_gs
38 : use m_ifc
39 : use m_symfind
40 :
41 : use defs_datatypes, only : pseudopotential_type
42 : use defs_abitypes, only : MPI_type
43 : use m_build_info, only : abinit_version
44 : use m_geometry, only : xcart2xred, xred2xcart
45 : use m_multibinit_dataset, only : multibinit_dtset_type
46 : use m_effective_potential, only : effective_potential_type
47 : use m_fit_polynomial_coeff, only : polynomial_coeff_writeXML, &
48 : fit_polynomial_coeff_fit, genereList, fit_polynomial_coeff_getPositive !,fit_polynomial_coeff_getCoeffBound
49 : use m_polynomial_coeff,only : polynomial_coeff_getNorder
50 : ! use m_pawang, only : pawang_type, pawang_free
51 : ! use m_pawrad, only : pawrad_type, pawrad_free
52 : ! use m_pawtab, only : pawtab_type, pawtab_nullify, pawtab_free
53 : ! use m_pawxmlps, only : paw_setup, ipsp2xml, rdpawpsxml, &
54 : !& paw_setup_copy, paw_setup_free, getecutfromxml
55 : use m_abihist
56 : use m_ewald
57 : use m_mpinfo, only : init_mpi_enreg,destroy_mpi_enreg
58 : use m_copy , only : alloc_copy
59 : use m_electronpositron, only : electronpositron_type
60 : use m_scfcv, only : scfcv_t, scfcv_run,scfcv_destroy
61 : use m_results_gs, only : results_gs_type
62 : use m_mover, only : mover
63 : use m_io_tools, only : get_unit, open_file
64 : use m_matrix, only: matr3inv,mati3inv
65 : use m_symtk, only: symatm
66 :
67 : implicit none
68 :
69 : private
70 : !!***
71 :
72 : public :: mover_effpot
73 : private :: checksymmetrygroup
74 : !!***
75 :
76 : contains
77 : !!***
78 :
79 : !!****f* ABINIT/mover_effpot
80 : !! NAME
81 : !! mover_effpot
82 : !!
83 : !! FUNCTION
84 : !! this routine is driver for using mover with effective potential
85 : !!
86 : !! INPUTS
87 : !! inp = input of multibinit
88 : !! effective_potential = effective potential of the reference structure
89 : !! option = flag for the option:
90 : !! -1 = Bound the anharmonic part
91 : !! 12 = NVT simulation
92 : !! 13 = NPT simulation
93 : !! OUTPUT
94 : !!
95 : !!
96 : !! SIDE EFFECTS
97 : !!
98 : !! NOTES
99 : !!
100 : !!
101 : !! SOURCE
102 :
103 44 : subroutine mover_effpot(inp,filnam,effective_potential,option,comm,hist)
104 :
105 : use m_scup_dataset, only : scup_dtset_type
106 :
107 : !Arguments --------------------------------
108 : !scalar
109 : integer, intent(in) :: option,comm
110 : !array
111 : type(multibinit_dtset_type),intent(inout) :: inp
112 : type(effective_potential_type),intent(inout) :: effective_potential
113 : character(len=fnlen),intent(in) :: filnam(15)
114 : type(abihist),optional,intent(inout):: hist
115 : !Local variables-------------------------------
116 : !scalar
117 : integer :: filetype,icoeff_bound,ii, unit_out,msym,nsym
118 : !integer :: iexit,initialized
119 : integer :: jj,kk,nproc,ncoeff,nmodels,ncoeff_bound,ncoeff_bound_tot,ncoeff_max
120 : integer :: model_bound,model_ncoeffbound,my_rank,ptgroupma,spgroup
121 : !integer :: mtypalch,,npsp,paw_size,type
122 : !integer,save :: paw_size_old=-1
123 : real(dp):: cutoff,freq_q,freq_b,qmass,bmass,tolsym
124 : ! real(dp):: time_q,time_b
125 : logical :: iam_master,isVused,isARused,readOnlyLast
126 : integer, parameter:: master=0
127 : logical :: verbose,writeHIST,file_opened
128 : !type
129 44 : type(scup_dtset_type) :: scup_inp
130 : !real(dp) :: cpui
131 : !character(len=8) :: codvsn
132 :
133 : !TEST_AM
134 : ! integer :: ia,mu,rand_seed = 5
135 : ! real(dp):: mass_ia,rescale_vel,sum_mass,v2gauss
136 : !TEST_AM
137 : ! Set array dimensions
138 : character(len=500) :: message
139 44 : type(MPI_type),target :: mpi_enreg
140 5632 : type(dataset_type),target :: dtset
141 44 : type(scfcv_t) :: scfcv_args
142 : type(datafiles_type),target :: dtfil
143 : integer,target :: zero_integer
144 44 : type(ab_xfh_type) :: ab_xfh
145 44 : type(results_gs_type),target :: results_gs
146 44 : type(pseudopotential_type),target :: psps
147 : !arrays
148 : !no_abirules
149 : integer :: sc_size(3),sc_size_TS(3)
150 44 : integer,pointer :: indsym(:,:,:)
151 44 : integer,allocatable :: listcoeff(:),listcoeff_bound(:,:),list_tmp(:),list_bound(:,:)
152 44 : integer,allocatable :: isPositive(:)
153 44 : integer,allocatable :: symrel(:,:,:),symrec(:,:,:)
154 : !integer,allocatable :: npwtot(:)
155 : real(dp) :: acell(3)
156 : !real(dp) :: ecut_tmp(3,2,10)
157 44 : real(dp),allocatable :: coeff_values(:,:)
158 44 : real(dp),pointer :: rhog(:,:),rhor(:,:)
159 44 : real(dp),allocatable :: tnons(:,:)
160 44 : real(dp),allocatable :: xred(:,:),xred_old(:,:),xcart(:,:)
161 44 : real(dp),allocatable :: gred(:,:),fcart(:,:)
162 44 : real(dp),allocatable :: vel(:,:)
163 : real(dp) :: vel_cell(3,3),rprimd(3,3)
164 : !real(dp) :: efield(3)
165 44 : type(polynomial_coeff_type),dimension(:),allocatable :: coeffs_all,coeffs_tmp,coeffs_bound
166 : character(len=fnlen) :: filename,md_hist_name
167 : character(len=fnlen) :: name_file
168 : !character(len=fnlen) :: filename_psp(3)
169 : type(electronpositron_type),pointer :: electronpositron
170 : ! type(pspheader_type),allocatable :: pspheads(:)
171 : ! type(pawrad_type),allocatable :: pawrad(:)
172 : ! type(pawtab_type),allocatable :: pawtab(:)
173 : ! type(args_gs_type) :: args_gs
174 : !type(wvl_data) :: wvl
175 : !type(pawang_type) :: pawang
176 : !type(scf_history_type) :: scf_history
177 44 : type(abihist) :: hist_tmp
178 :
179 :
180 : !******************************************************************
181 :
182 : !MPI variables
183 44 : nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
184 44 : iam_master = (my_rank == master)
185 :
186 44 : write(message, '(a,(80a),a)') ch10,&
187 3608 : & ('=',ii=1,80),ch10
188 44 : call wrtout(ab_out,message,'COLL')
189 44 : call wrtout(std_out,message,'COLL')
190 :
191 : !*******************************************************************
192 : ! 1 Generate and check supercell for the dynamics
193 : !*******************************************************************
194 :
195 : !a new supercell is compute
196 : !Initialisaton of variable
197 44 : if(option == -1.or.option == -2) then
198 : ! Bound process option
199 0 : sc_size(:) = inp%bound_cell
200 44 : else if(option == -3) then
201 : ! Heff option
202 0 : sc_size(:) = (/1,1,1/)
203 : else
204 : ! Normal dynamics
205 176 : sc_size(:) = inp%ncell
206 : end if
207 :
208 44 : if(option/=0)then
209 :
210 : ! acell is always set to one, only rprimd is used for the effective potential
211 176 : acell = one
212 572 : rprimd = effective_potential%crystal%rprimd
213 :
214 132 : ABI_MALLOC(xred,(3,effective_potential%crystal%natom))
215 88 : ABI_MALLOC(xcart,(3,effective_potential%crystal%natom))
216 :
217 : ! convert new xcart
218 : call xcart2xred(effective_potential%crystal%natom,effective_potential%crystal%rprimd,&
219 44 : & effective_potential%crystal%xcart,xred)
220 44 : call xred2xcart(effective_potential%crystal%natom, rprimd, xcart, xred)
221 : ! Generate supercell for the simulation
222 44 : call effective_potential_setSupercell(effective_potential,comm,ncell=sc_size)
223 :
224 44 : ABI_FREE(xred)
225 44 : ABI_FREE(xcart)
226 :
227 : !***************************************************************
228 : !1 Convert some parameters into the structures used by mover.F90
229 : !***************************************************************
230 : !NOTE:ARGUMENTS OF MOVER SHOULD BE CLEAN
231 : ! We may just need to provide AB_MOVER wich is the main object
232 : ! for mover and set scfcv_args as an optional and depending on
233 : ! the kind of calculation (abinit or multibinit), we provide
234 : !! to mover scfcv_args or effective_potential...
235 : !***************************************************************
236 : ! Free dtset
237 44 : call dtset%free()
238 :
239 : ! Set mpi_eng
240 44 : mpi_enreg%comm_cell = comm
241 44 : mpi_enreg%me = my_rank
242 :
243 : ! Set the abinit dataset for mover with fake values
244 : ! Scalar
245 44 : dtset%dmft_entropy = 0
246 44 : dtset%nctime = inp%nctime ! NetCdf TIME between output of molecular dynamics information
247 44 : dtset%delayperm = 0 ! DELAY between trials to PERMUTE atoms
248 44 : dtset%dilatmx = 1.5 ! DILATation : MaXimal value
249 44 : dtset%chkdilatmx = 0 ! No check on dilatmx is needed in multibilint
250 44 : dtset%diismemory = 8 ! Direct Inversion in the Iterative Subspace MEMORY
251 44 : dtset%friction = 0.0001d0 ! internal FRICTION coefficient
252 44 : dtset%goprecon = 0 ! Geometry Optimization PREconditioner equations
253 44 : dtset%istatr = 0 ! Integer for STATus file SHiFT
254 44 : dtset%jellslab = 0 ! include a JELLium SLAB in the cell
255 44 : dtset%mqgrid = 0 ! Maximum number of Q-space GRID points for pseudopotentials
256 44 : dtset%mqgriddg = 0 ! Maximum number of Q-wavevectors for the 1-dimensional GRID
257 : ! for the Double Grid in PAW
258 44 : dtset%mdwall = 10000d0 ! Molecular Dynamics WALL location
259 44 : dtset%ntypalch = 0 ! Number of TYPe of atoms that are "ALCHemical"
260 44 : dtset%natom = effective_potential%supercell%natom
261 44 : dtset%ntypat = effective_potential%crystal%ntypat
262 44 : dtset%npspalch = effective_potential%crystal%ntypat
263 44 : dtset%nconeq = 0 ! Number of CONstraint EQuations
264 44 : dtset%noseinert = 1.d-5 ! NOSE INERTia factor
265 44 : dtset%nnos = inp%nnos ! Number of nose masses Characteristic
266 176 : do ii=1,3 ! Only copy the diagonal part
267 : ! might be adapt with ngqpt(9) instead of ngqpt(3) because is needed in wght9.f
268 176 : dtset%ph_ngqpt(ii) = inp%ngqpt((ii))
269 : end do
270 44 : dtset%ph_nqshift = inp%nqshft
271 44 : dtset%prtxml = 0 ! print the xml
272 44 : dtset%signperm = 1 ! SIGN of PERMutation potential
273 44 : dtset%strprecon = inp%strprecon ! STRess PRECONditioner
274 176 : dtset%supercell_latt(:) = 0
275 176 : do ii=1,3
276 176 : dtset%supercell_latt(ii) = sc_size(ii)
277 : end do
278 44 : dtset%tolmxf = inp%tolmxf
279 44 : dtset%tsmear = 0.009500446 !
280 44 : dtset%vis = 100 ! VIScosity
281 44 : dtset%usewvl = 0 !
282 44 : dtset%useylm = 0 !
283 :
284 : ! if(option == -3) then
285 : ! write(message,'(a)')' Read the DDB file to fill the dtset array'
286 : ! call wrtout(std_out,message,"COLL")
287 : ! ! Copy real informtions from the ddb
288 : ! call effective_potential_file_getType(filnam(3),type)
289 : ! if(type /= 1) then
290 : ! write(message, '(5a)' )&
291 : ! & ' You need to provide DDB file in the input to compute ahnarmonic',ch10,&
292 : ! & ' part of effective Hamiltionian',ch10,&
293 : ! & 'Action: add DDB file in the inputs'
294 : ! ABI_BUG(message)
295 : ! end if
296 : ! call ddb_to_dtset(comm, dtset,filnam(3),psps)
297 : ! ABI_MALLOC(dtset%kptns,(3,dtset%nkpt))
298 : ! dtset%kptns(:,:) = dtset%kpt(:,:)
299 : ! ABI_MALLOC(dtset%istwfk,(dtset%nkpt))
300 : ! dtset%istwfk(:) = 1
301 :
302 : ! else
303 : ! Need to init some values
304 : ! dtset%nsym = 1 ! Number of SYMmetry operations
305 : ! ABI_MALLOC(symrel,(3,3,dtset%nsym))
306 : ! symrel = reshape((/1,0,0,0,1,0,0,0,1/),shape(symrel))
307 : ! call alloc_copy(symrel,dtset%symrel)
308 : ! ABI_MALLOC(tnons,(3,dtset%nsym))
309 : ! tnons = zero
310 : ! call alloc_copy(tnons,dtset%tnons)
311 44 : call alloc_copy(effective_potential%supercell%typat,dtset%typat)
312 44 : call alloc_copy(effective_potential%crystal%znucl,dtset%znucl)
313 : ! ABI_FREE(symrel)
314 : ! ABI_FREE(tnons)
315 : ! end if
316 : !Find symmetry for simulation
317 44 : if(inp%dyn_chksym == 1)then
318 : write(message,'(2a)')&
319 1 : & ' Check Symmetry of Start Structure and Impose it for dynamics run',ch10
320 1 : call wrtout(std_out,message,"COLL")
321 1 : call wrtout(ab_out,message,"COLL")
322 : !read md_hist
323 1 : if(inp%restartxf < 0)then
324 1 : md_hist_name=trim(filnam(2))//'_HIST.nc'
325 : write(message,'(3a)')&
326 1 : & ' Restart from external structure stored in file: ',trim(md_hist_name),ch10
327 1 : call wrtout(std_out,message,"COLL")
328 1 : call wrtout(ab_out,message,"COLL")
329 1 : isVused = .true.
330 1 : isARused = .true.
331 1 : readOnlyLast = .true.
332 1 : call read_md_hist(md_hist_name,hist_tmp,isVused,isARused,readOnlyLast)
333 :
334 :
335 :
336 : write(message,'(2a)')&
337 1 : & ' Map external structure to internal ordering of reference structure: ',ch10
338 1 : call wrtout(std_out,message,"COLL")
339 1 : call effective_potential_file_mapHistToRef(effective_potential,hist_tmp,comm,verbose=.True.) ! Map Hist to Ref to order atoms
340 : endif
341 4 : msym = 96*PRODUCT(sc_size)
342 3 : ABI_MALLOC(symrel,(3,3,msym))
343 2 : ABI_MALLOC(symrec,(3,3,msym))
344 3 : ABI_MALLOC(tnons,(3,msym))
345 1 : tolsym = inp%dyn_tolsym
346 : call checksymmetrygroup(hist_tmp%rprimd,hist_tmp%xred,&
347 : & effective_potential%supercell%typat,msym,effective_potential%supercell%natom,&
348 1 : & ptgroupma,spgroup,symrel,tnons,nsym,tolsym)
349 1 : dtset%nsym = nsym
350 17 : do ii=1,nsym
351 17 : call mati3inv(symrel(:,:,ii),symrec(:,:,ii))
352 : end do
353 : !Get Indsym
354 4 : ABI_MALLOC(indsym,(4,dtset%nsym,dtset%natom))
355 1 : call symatm(indsym,dtset%natom,dtset%nsym,symrec,tnons,tolsym,effective_potential%supercell%typat,hist_tmp%xred)
356 1 : call alloc_copy(symrel,dtset%symrel)
357 1 : call alloc_copy(tnons,dtset%tnons)
358 1 : scfcv_args%indsym => indsym
359 1 : call abihist_free(hist_tmp)
360 2 : ABI_FREE(symrec)
361 : else
362 43 : dtset%nsym = 1 ! Number of SYMmetry operations
363 43 : ABI_MALLOC(symrel,(3,3,dtset%nsym))
364 172 : symrel = reshape((/1,0,0,0,1,0,0,0,1/),shape(symrel))
365 43 : call alloc_copy(symrel,dtset%symrel)
366 129 : ABI_MALLOC(tnons,(3,dtset%nsym))
367 215 : tnons = zero
368 43 : call alloc_copy(tnons,dtset%tnons)
369 172 : ABI_MALLOC(indsym,(4,dtset%nsym,dtset%natom))
370 28771 : indsym = 0
371 43 : scfcv_args%indsym => indsym
372 : endif
373 :
374 44 : ABI_FREE(symrel)
375 44 : ABI_FREE(tnons)
376 :
377 : !array
378 132 : ABI_MALLOC(dtset%iatfix,(3,dtset%natom)) ! Indices of AToms that are FIXed
379 19320 : dtset%iatfix = inp%iatfix
380 176 : dtset%goprecprm(:) = zero !Geometry Optimization PREconditioner PaRaMeters equations
381 132 : ABI_MALLOC(dtset%prtatlist,(dtset%natom)) !PRinT by ATom LIST of ATom
382 4852 : dtset%prtatlist(:) = 0
383 176 : ABI_MALLOC(dtset%mixalch_orig,(dtset%npspalch,dtset%ntypalch,1))
384 88 : dtset%mixalch_orig(:,:,:)=zero
385 132 : ABI_MALLOC(dtset%ph_qshift,(3,dtset%ph_nqshift))
386 264 : dtset%ph_qshift = inp%q1shft
387 44 : dtset%hmctt = inp%hmctt
388 44 : dtset%hmcsst = inp%hmcsst
389 44 : if(option > 0)then
390 44 : verbose = .TRUE.
391 44 : writeHIST = .TRUE.
392 44 : dtset%dtion = inp%dtion ! Delta Time for IONs
393 44 : dtset%ionmov = inp%dynamics ! Number for the dynamic
394 44 : dtset%ntime = inp%ntime ! Number of TIME steps
395 44 : dtset%optcell = inp%optcell ! OPTimize the CELL shape and dimensions Characteristic
396 44 : dtset%restartxf = inp%restartxf ! RESTART from (X,F) history
397 44 : dtset%mdtemp(1) = inp%temperature !Molecular Dynamics Temperatures
398 44 : dtset%mdtemp(2) = inp%temperature !Molecular Dynamics Temperatures
399 44 : dtset%strfact = inp%strfact ! STRess FACTor
400 308 : dtset%strtarget(1:6) = -1 * inp%strtarget(1:6) / 29421.033d0 ! STRess TARGET
401 0 : else if(option == -1.or.option == -2) then
402 : ! Set default for the fit
403 0 : verbose = .false.
404 0 : writeHIST = .false.
405 0 : dtset%restartxf = 0 ! RESTART from (X,F) history
406 0 : dtset%dtion = inp%dtion ! Delta Time for IONs
407 0 : dtset%ionmov = 13 ! Number for the dynamic
408 0 : dtset%ntime = inp%bound_step ! Number of TIME steps
409 0 : dtset%optcell = 2 ! OPTimize the CELL shape and dimensions Characteristic
410 0 : dtset%mdtemp(1) = inp%bound_temp !Molecular Dynamics Temperatures
411 0 : dtset%mdtemp(2) = inp%bound_temp !Molecular Dynamics Temperatures
412 0 : dtset%strfact = 100.0d0
413 0 : dtset%strtarget(1:6) = zero
414 : end if
415 :
416 : ! Set the barostat and thermonstat if ionmov == 13
417 44 : if(dtset%ionmov == 13 .or. dtset%ionmov == 25)then
418 :
419 : ! Select frequency of the barostat as a function of temperature
420 : ! For small temperature, we need huge barostat and inversely
421 : ! if(dtset%mdtemp(1) <= 10) then
422 : ! freq_q = 0.002
423 : ! freq_b = 0.0002
424 : ! else if(dtset%mdtemp(1) <= 50) then
425 : ! freq_q = 0.02
426 : ! freq_b = 0.002
427 : ! else if(dtset%mdtemp(1) > 50.and.dtset%mdtemp(1) < 300) then
428 : ! freq_q = 0.1
429 : ! freq_b = 0.01
430 : ! else
431 : ! freq_q = 0.2
432 : ! freq_b = 0.02
433 : ! end if
434 :
435 : !TEST_AM_old way
436 : ! freq_q = 0.1
437 : ! freq_b = 0.01
438 : ! qmass = dtset%natom* kb_THzK * dtset%mdtemp(1) / (freq_q**2)
439 : ! bmass = dtset%natom* kb_THzK * dtset%mdtemp(1) / (freq_b**2)
440 : !TEST_AM
441 :
442 :
443 : !TEST_AM
444 26 : freq_q = 800 / Ha_cmm1
445 26 : freq_b = 800 / Ha_cmm1
446 26 : qmass = 10 * dtset%natom * kb_HaK * dtset%mdtemp(1) / (freq_q**2)
447 26 : bmass = 10000*qmass
448 : !TEST_AM
449 :
450 :
451 26 : if(dtset%nnos==0) then
452 26 : dtset%nnos = 1
453 26 : ABI_MALLOC(dtset%qmass,(dtset%nnos))
454 52 : dtset%qmass(:) = qmass
455 : write(message,'(3a,F30.10,a)')&
456 26 : & ' WARNING: nnos is set to zero in the input',ch10,&
457 52 : & ' value by default for qmass: ',dtset%qmass(:),ch10
458 26 : if(verbose)call wrtout(std_out,message,"COLL")
459 : else
460 0 : ABI_MALLOC(dtset%qmass,(dtset%nnos)) ! Q thermostat mass
461 0 : dtset%qmass(:) = inp%qmass(:)
462 : end if
463 26 : if (abs(inp%bmass) < tol10) then
464 26 : dtset%bmass = bmass
465 : write(message,'(3a,F30.10,a)')&
466 26 : & ' WARNING: bmass is set to zero in the input',ch10,&
467 52 : & ' value by default for bmass: ',dtset%bmass,ch10
468 26 : if(verbose)call wrtout(std_out,message,"COLL")
469 : else
470 0 : dtset%bmass = inp%bmass ! Barostat mass
471 : end if
472 : end if
473 :
474 : ! Do a check
475 44 : if(dtset%ionmov == 27)then
476 0 : call effective_potential_file_getType(filnam(3),filetype)
477 0 : if(filetype /= 1)then
478 : write(message, '(5a)' )&
479 0 : & ' The file ',trim(filnam(3)),' is not a DDB',ch10,&
480 0 : & ' It is not compatible with ionmov 27'
481 0 : ABI_ERROR(message)
482 : end if
483 :
484 : end if
485 :
486 : ! set psps
487 44 : psps%useylm = dtset%useylm
488 :
489 : ! initialisation of results_gs
490 44 : call results_gs%init(dtset%natom,1,1)
491 :
492 : ! Set the pointers of scfcv_args
493 44 : zero_integer = 0
494 44 : scfcv_args%dtset => dtset
495 : !ABI_MALLOC(indsym,(4,dtset%nsym,dtset%natom))
496 : !indsym = 0
497 : !scfcv_args%indsym => indsym
498 44 : scfcv_args%mpi_enreg => mpi_enreg
499 44 : scfcv_args%ndtpawuj => zero_integer
500 44 : scfcv_args%results_gs => results_gs
501 44 : scfcv_args%psps => psps
502 : ! Set other arguments of the mover.F90 routines
503 :
504 : ! Set the dffil structure
505 132 : dtfil%filnam_ds(1:2)=filnam(1:2)
506 44 : dtfil%filnam_ds(3)=""
507 44 : dtfil%filnam_ds(4)=filnam(2)
508 44 : dtfil%filstat='_STATUS'
509 44 : nullify (electronpositron)
510 44 : ABI_MALLOC(rhog,(2,1))
511 44 : ABI_MALLOC(rhor,(2,1))
512 :
513 : ! Initialize xf history (should be put in inwffil)
514 : ! Not yet implemented for ionmov 2 3 10 11 22 (memory problem...)
515 : ! ab_xfh%mxfh=(ab_xfh%nxfh-dtset%restartxf+1)+dtset%ntime+5
516 44 : ab_xfh%nxfh = 0
517 44 : ab_xfh%mxfh = 1
518 176 : ABI_MALLOC(ab_xfh%xfhist,(3,dtset%natom+4,2,ab_xfh%mxfh))
519 44 : if (any((/3,10,11/)==dtset%ionmov)) then
520 : write(message, '(3a)' )&
521 0 : & ' This dynamics can not be used with effective potential',ch10,&
522 0 : & 'Action: correct dynamics input'
523 0 : ABI_BUG(message)
524 : end if
525 :
526 : !Get SCALE-UP INPUT
527 :
528 44 : scup_inp = inp%scup_dtset
529 :
530 :
531 : !***************************************************************
532 : !2 initialization of the structure for the dynamics
533 : !***************************************************************
534 :
535 44 : if (allocated(dtset%rprimd_orig)) then
536 0 : ABI_FREE(dtset%rprimd_orig)
537 : end if
538 44 : ABI_MALLOC(dtset%rprimd_orig,(3,3,1))
539 572 : dtset%rprimd_orig(:,:,1) = effective_potential%supercell%rprimd
540 :
541 :
542 132 : ABI_MALLOC(xred,(3,dtset%natom))
543 88 : ABI_MALLOC(xred_old,(3,dtset%natom))
544 88 : ABI_MALLOC(vel,(3,dtset%natom))
545 88 : ABI_MALLOC(gred,(3,dtset%natom))
546 88 : ABI_MALLOC(fcart,(3,dtset%natom))
547 :
548 : call xcart2xred(dtset%natom,effective_potential%supercell%rprimd,&
549 44 : & effective_potential%supercell%xcart,xred)
550 :
551 19320 : xred_old = xred
552 44 : vel_cell(:,:) = zero
553 19276 : vel(:,:) = zero
554 :
555 : !*********************************************************
556 : !4 Call main routine for the bound process,
557 : ! monte carlo / molecular dynamics / project
558 : !*********************************************************
559 44 : if(option > 0)then
560 : !*************************************************************
561 : ! call mover in case of NPT or NVT simulation
562 : !*************************************************************
563 3608 : write(message, '((80a),3a)' ) ('-',ii=1,80), ch10,&
564 3652 : & '-Monte Carlo / Molecular Dynamics ',ch10
565 :
566 : ! Marcus: if wanted analyze anharmonic terms of effective potential &&
567 : ! and print anharmonic contribution to file anharmonic_energy_terms.out
568 : ! Open File and write header
569 44 : ncoeff = effective_potential%anharmonics_terms%ncoeff
570 44 : name_file='MD'
571 44 : if(inp%analyze_anh_pot == 1)then
572 : call effective_potential_writeAnhHead(ncoeff,name_file,&
573 0 : & effective_potential%anharmonics_terms)
574 : end if
575 :
576 44 : call wrtout(ab_out,message,'COLL')
577 44 : call wrtout(std_out,message,'COLL')
578 : call mover(scfcv_args,ab_xfh,acell,effective_potential%crystal%amu,dtfil,electronpositron,&
579 : & rhog,rhor,dtset%rprimd_orig,vel,vel_cell,xred,xred_old,&
580 : & effective_potential=effective_potential,filename_ddb=filnam(3),&
581 : & verbose=verbose,verbose_cycle=dtset%ionmov/=25,writeHIST=writeHIST,scup_dtset=scup_inp,sc_size=sc_size(:),&
582 44 : & multibinit_dtset=inp)
583 44 : INQUIRE(FILE='MD_anharmonic_terms_energy.dat',OPENED=file_opened,number=unit_out)
584 44 : if(file_opened) close(unit_out)
585 0 : else if(option== -1.or.option==-2)then
586 : !*************************************************************
587 : ! Try to bound the model
588 : !*************************************************************
589 0 : write(message, '((80a),4a)' ) ('-',ii=1,80), ch10,&
590 0 : & ' Try to bound the model',ch10,' Check if the model is bounded or not'
591 0 : call wrtout(ab_out,message,'COLL')
592 0 : call wrtout(std_out,message,'COLL')
593 :
594 : ! Try the model
595 : call mover(scfcv_args,ab_xfh,acell,effective_potential%crystal%amu,dtfil,electronpositron,&
596 : & rhog,rhor,dtset%rprimd_orig,vel,vel_cell,xred,xred_old,&
597 : & effective_potential=effective_potential,verbose=verbose,verbose_cycle=dtset%ionmov/=25,writeHIST=writeHIST,&
598 0 : & multibinit_dtset=inp )
599 :
600 0 : write(message, '(a)' ) ' => The model'
601 0 : if(effective_potential%anharmonics_terms%bounded)then
602 0 : write(message, '(2a)' ) trim(message),' is bound'
603 0 : call wrtout(std_out,message,'COLL')
604 0 : call wrtout(ab_out,message,'COLL')
605 : else
606 0 : write(message, '(2a)' ) trim(message),' is not bound'
607 0 : call wrtout(std_out,message,'COLL')
608 :
609 :
610 0 : if(option==-2)then
611 :
612 : ! Fill the list for the fixcoeff input of the fit_polynomial_coeff_fit routine
613 : ! Store the number of coefficients before adding other coeff for the bounding
614 0 : ncoeff = effective_potential%anharmonics_terms%ncoeff
615 0 : ABI_MALLOC(listcoeff,(ncoeff))
616 0 : do ii=1,ncoeff
617 0 : listcoeff(ii)=ii
618 : end do
619 :
620 0 : write(message, '(2a)')ch10,' Generate the list of addionnal terms with the fit process...'
621 0 : call wrtout(std_out,message,'COLL')
622 :
623 : ! Get the additional coeff
624 : call fit_polynomial_coeff_fit(effective_potential,(/0/),listcoeff,hist,1,&
625 : & inp%bound_rangePower,0,inp%bound_maxCoeff,ncoeff,inp%fit_nimposecoeff,inp%fit_imposecoeff,&
626 : & 1,comm,cutoff_in=inp%bound_cutoff,&
627 : & max_power_strain=2,verbose=.true.,positive=.true.,spcoupling=inp%bound_SPCoupling==1,&
628 : & anharmstr=inp%bound_anhaStrain==1,only_even_power=.true.,fit_on=inp%fit_on,sel_on=inp%sel_on, &
629 0 : & max_nbody=inp%fit_max_nbody, drop_rate=0.0_dp, ncoeff_per_cycle=1, fit_weight_T=inp%fit_weight_T)
630 :
631 : ! Store the max number of coefficients after the fit process
632 0 : ncoeff_max = effective_potential%anharmonics_terms%ncoeff
633 : ! Store all the coefficients in coeffs_all
634 0 : ABI_MALLOC(coeffs_all,(ncoeff_max))
635 0 : ABI_MALLOC(coeffs_tmp,(ncoeff_max))
636 0 : do ii=1,ncoeff_max
637 : call polynomial_coeff_init(&
638 : & effective_potential%anharmonics_terms%coefficients(ii)%coefficient,&
639 : & effective_potential%anharmonics_terms%coefficients(ii)%nterm,&
640 : & coeffs_all(ii),&
641 : & effective_potential%anharmonics_terms%coefficients(ii)%terms,&
642 : & effective_potential%anharmonics_terms%coefficients(ii)%name,&
643 0 : & check=.false.)
644 : end do
645 :
646 0 : model_ncoeffbound = 0
647 0 : do ii=1,ncoeff_max-ncoeff
648 :
649 0 : write(message, '(5a,I0,a)')ch10,'--',ch10,' Try to bound the model ',&
650 0 : & 'with ', ii,' additional term'
651 0 : if(ii>1)write(message,'(2a)') trim(message),'s'
652 0 : call wrtout(std_out,message,'COLL')
653 0 : call wrtout(ab_out,message,'COLL')
654 :
655 : ! Copy the new model in coeffs_tmp(jj)
656 : ! Free the coeffs_tmp array before
657 0 : do jj=1,ncoeff_max
658 0 : call polynomial_coeff_free(coeffs_tmp(jj))
659 : end do
660 :
661 0 : do jj=1,ncoeff+ii
662 : call polynomial_coeff_init(&
663 : & coeffs_all(jj)%coefficient,&
664 : & coeffs_all(jj)%nterm,&
665 : & coeffs_tmp(jj),&
666 : & coeffs_all(jj)%terms,&
667 : & coeffs_all(jj)%name,&
668 0 : & check=.false.)
669 : end do
670 :
671 0 : model_ncoeffbound = ii
672 :
673 : ! Reset the simulation and set the coefficients of the model
674 0 : call effective_potential_setCoeffs(coeffs_tmp(1:ncoeff+ii),effective_potential,ncoeff+ii)
675 : call fit_polynomial_coeff_fit(effective_potential,(/0/),(/0/),hist,0,(/0,0/),0,0,&
676 : & -1,inp%fit_nimposecoeff,inp%fit_imposecoeff,1,comm,verbose=.true.,positive=.false.,&
677 : & fit_on=inp%fit_on,sel_on=inp%sel_on, max_nbody=inp%fit_max_nbody, drop_rate=0.0_dp, &
678 0 : & ncoeff_per_cycle=1, fit_weight_T=inp%fit_weight_T)
679 0 : call effective_potential_setSupercell(effective_potential,comm,ncell=sc_size)
680 0 : dtset%rprimd_orig(:,:,1) = effective_potential%supercell%rprimd
681 0 : acell(:) = one
682 : call xcart2xred(dtset%natom,effective_potential%supercell%rprimd,&
683 0 : & effective_potential%supercell%xcart,xred)
684 0 : xred_old = xred
685 0 : vel_cell(:,:) = zero
686 0 : vel(:,:) = zero
687 0 : gred(:,:) = zero
688 0 : fcart(:,:) = zero
689 :
690 : ! Run mover to check if the model is bound
691 : call mover(scfcv_args,ab_xfh,acell,effective_potential%crystal%amu,dtfil,electronpositron,&
692 : & rhog,rhor,dtset%rprimd_orig,vel,vel_cell,xred,xred_old,&
693 0 : & effective_potential=effective_potential,verbose=verbose,verbose_cycle=dtset%ionmov/=25,writeHIST=writeHIST)
694 0 : if(.not.effective_potential%anharmonics_terms%bounded)then
695 0 : write(message, '(2a)' ) ' => The model is not bounded'
696 : else
697 0 : write(message, '(2a)' ) ' => The model is bounded'
698 : end if
699 0 : call wrtout(std_out,message,'COLL')
700 :
701 : ! Exit if the model is bounded
702 0 : if(effective_potential%anharmonics_terms%bounded) exit
703 :
704 : end do
705 :
706 : else
707 : ! Get the list of possible coefficients to bound the model
708 0 : cutoff = zero
709 0 : do ii=1,3
710 0 : cutoff = cutoff + effective_potential%crystal%rprimd(ii,ii)
711 : end do
712 0 : cutoff = cutoff / 3.0
713 :
714 : ! call fit_polynomial_coeff_getCoeffBound(effective_potential,coeffs_bound,&
715 : !& hist,ncoeff_bound,comm,verbose=.true.)
716 :
717 : !TEST_AM!
718 0 : sc_size_TS = (/2,2,2/)
719 : call polynomial_coeff_getNorder(coeffs_bound,effective_potential%crystal,cutoff,&
720 : & ncoeff_bound,ncoeff_bound_tot,inp%bound_rangePower,inp%bound_rangePower(2),2,sc_size_TS,&
721 : & comm,anharmstr=inp%bound_anhaStrain==1,&
722 : & spcoupling=inp%bound_SPCoupling==1,verbose=.false.,distributed=.false.,&
723 0 : & only_even_power=.true.,only_odd_power=.false., max_nbody=inp%fit_max_nbody)
724 :
725 0 : if(iam_master)then
726 0 : filename=trim(filnam(2))//"_boundcoeff.xml"
727 0 : call polynomial_coeff_writeXML(coeffs_bound,ncoeff_bound,filename=filename,newfile=.true.)
728 : end if
729 : ! wait
730 0 : call xmpi_barrier(comm)
731 : ! Store all the initial coefficients
732 0 : ncoeff = effective_potential%anharmonics_terms%ncoeff
733 0 : ABI_MALLOC(coeffs_all,(ncoeff+ncoeff_bound))
734 0 : do ii=1,ncoeff
735 : call polynomial_coeff_init(effective_potential%anharmonics_terms%coefficients(ii)%coefficient,&
736 : & effective_potential%anharmonics_terms%coefficients(ii)%nterm,&
737 : & coeffs_all(ii),&
738 : & effective_potential%anharmonics_terms%coefficients(ii)%terms,&
739 : & effective_potential%anharmonics_terms%coefficients(ii)%name,&
740 0 : & check=.false.)
741 : end do
742 :
743 0 : do ii=1,ncoeff_bound
744 : call polynomial_coeff_init(coeffs_bound(ii)%coefficient,&
745 : & coeffs_bound(ii)%nterm,&
746 : & coeffs_all(ncoeff+ii),&
747 : & coeffs_bound(ii)%terms,&
748 : & coeffs_bound(ii)%name,&
749 0 : & check=.false.)
750 : end do
751 :
752 : ! Copy the fixed coefficients from the model (without bound coeff)
753 0 : ncoeff = effective_potential%anharmonics_terms%ncoeff
754 0 : ABI_MALLOC(coeffs_tmp,(ncoeff+ncoeff_bound))
755 0 : do ii=1,ncoeff
756 : call polynomial_coeff_init(effective_potential%anharmonics_terms%coefficients(ii)%coefficient,&
757 : & effective_potential%anharmonics_terms%coefficients(ii)%nterm,&
758 : & coeffs_tmp(ii),&
759 : & effective_potential%anharmonics_terms%coefficients(ii)%terms,&
760 : & effective_potential%anharmonics_terms%coefficients(ii)%name,&
761 0 : & check=.false.)
762 : end do
763 :
764 0 : ncoeff_max = ncoeff+ncoeff_bound
765 0 : ABI_MALLOC(listcoeff,(ncoeff_max))
766 0 : listcoeff = 0
767 0 : do jj=1,ncoeff
768 0 : listcoeff(jj) = jj
769 : end do
770 :
771 0 : model_bound = 0
772 0 : model_ncoeffbound = 0
773 :
774 0 : do ii=2,inp%bound_maxCoeff
775 : ! Compute the number of possible combination
776 0 : nmodels = 1
777 0 : ABI_MALLOC(list_bound,(nmodels,ii))
778 0 : ABI_MALLOC(list_tmp,(ii))
779 0 : list_bound = 0; list_tmp = 0; kk = 0; jj = 1
780 :
781 : ! Generate the list of possible combinaison 1st count
782 0 : call genereList(kk,jj,ii,ncoeff_bound,list_tmp,list_bound,nmodels,.false.)
783 0 : nmodels = kk
784 :
785 0 : write(message, '(5a,I0,a,I0,a)')ch10,'--',ch10,' Try to bound the model ',&
786 0 : & 'with ', ii,' additional positive terms (',nmodels,') possibilities'
787 0 : call wrtout(std_out,message,'COLL')
788 :
789 : ! allocate and generate combinaisons
790 0 : ABI_FREE(list_bound)
791 0 : ABI_FREE(list_tmp)
792 0 : ABI_MALLOC(coeff_values,(nmodels,ncoeff+ii))
793 0 : ABI_MALLOC(listcoeff_bound,(nmodels,ncoeff+ii))
794 0 : ABI_MALLOC(list_bound,(nmodels,ii))
795 0 : ABI_MALLOC(list_tmp,(ii))
796 0 : ABI_MALLOC(isPositive,(nmodels))
797 0 : list_bound = 0; listcoeff_bound = 0; list_tmp = 0; isPositive = 0; kk = 0; jj = 1
798 0 : call genereList(kk,jj,ii,ncoeff_bound,list_tmp,list_bound,nmodels,.true.)
799 : ! Generate the models
800 0 : do jj=1,nmodels
801 0 : listcoeff_bound(jj,1:ncoeff) = listcoeff(1:ncoeff)
802 0 : listcoeff_bound(jj,ncoeff+1:ncoeff+ii) = list_bound(jj,:) + ncoeff
803 : end do
804 :
805 : ! Reset the simulation
806 0 : call effective_potential_setCoeffs(coeffs_all,effective_potential,ncoeff+ncoeff_bound)
807 : call fit_polynomial_coeff_getPositive(effective_potential,hist,coeff_values,&
808 0 : & isPositive,listcoeff_bound,ncoeff+ii,ncoeff,nmodels,comm,verbose=.false., fit_weight_T=inp%fit_weight_T)
809 0 : if(all(isPositive == 0)) then
810 0 : write(message, '(5a,I0,a)')ch10,'--',ch10,' No possible model ',&
811 0 : & 'with ', ii,' additional terms found'
812 0 : call wrtout(std_out,message,'COLL')
813 : else
814 :
815 0 : do jj=1,nmodels
816 0 : if(isPositive(jj) == 1 .and. all(abs(coeff_values(jj,:)) < 1.0E5)) then
817 0 : write(message, '(2a,I0,a)') ch10,' The model number ',jj,' ['
818 0 : do kk=1,ncoeff+ii
819 0 : if(kk<ncoeff+ii)then
820 0 : write(message, '(a,I0,a)') trim(message),listcoeff_bound(jj,kk),','
821 : else
822 0 : write(message, '(a,I0)') trim(message),listcoeff_bound(jj,kk)
823 : end if
824 : end do
825 0 : write(message, '(2a)') trim(message),'] is positive'
826 0 : call wrtout(std_out,message,'COLL')
827 0 : write(message, '(2a,I0,a)') ' Check if the model ',&
828 0 : & 'number ', jj,' is bounded...'
829 0 : call wrtout(std_out,message,'COLL')
830 :
831 : ! Set the coefficients of the model
832 0 : do kk=1,ncoeff+ii
833 0 : if(kk<=ncoeff)then
834 : ! just set the values of the coefficient
835 0 : call polynomial_coeff_setCoefficient(coeff_values(jj,kk),coeffs_tmp(kk))
836 0 : write(message, '(a,I0,a,ES19.10,2a)') ' Set the value of the coefficient ',kk,&
837 0 : & ' =>',coeff_values(jj,kk),' ',trim(coeffs_tmp(kk)%name)
838 0 : call wrtout(std_out,message,'COLL')
839 :
840 : else
841 : ! Set the good coefficient
842 0 : icoeff_bound = listcoeff_bound(jj,kk)-ncoeff ! need to remove ncoeff value
843 : write(message, '(a,I0,a,I0,a,ES19.10,2a)')&
844 0 : & ' Set the value of the coefficient ',kk,' (',icoeff_bound,&
845 0 : & ') =>',coeff_values(jj,kk),&
846 0 : & ' ',trim(coeffs_bound(icoeff_bound)%name)
847 0 : call wrtout(std_out,message,'COLL')
848 0 : call polynomial_coeff_free(coeffs_tmp(kk))
849 : call polynomial_coeff_init(coeff_values(jj,kk),&
850 : & coeffs_bound(icoeff_bound)%nterm,&
851 : & coeffs_tmp(kk),&
852 : & coeffs_bound(icoeff_bound)%terms,&
853 : & coeffs_bound(icoeff_bound)%name,&
854 0 : & check=.false.)
855 :
856 : end if
857 : end do
858 :
859 : ! Reset the simulation and set the coefficients of the model
860 : call effective_potential_setCoeffs(coeffs_tmp(1:ncoeff+ii),effective_potential,&
861 0 : & ncoeff+ii)
862 : call fit_polynomial_coeff_fit(effective_potential,(/0/),(/0/),hist,0,(/0,0/),1,0,&
863 : & -1,inp%fit_nimposecoeff,inp%fit_imposecoeff,1,comm,verbose=.false.,positive=.false., &
864 : & max_nbody=inp%fit_max_nbody, drop_rate=0.0_dp, ncoeff_per_cycle=1, &
865 0 : & fit_weight_T=inp%fit_weight_T)
866 0 : call effective_potential_setSupercell(effective_potential,comm,ncell=sc_size)
867 0 : dtset%rprimd_orig(:,:,1) = effective_potential%supercell%rprimd
868 0 : acell(:) = one
869 : call xcart2xred(dtset%natom,effective_potential%supercell%rprimd,&
870 0 : & effective_potential%supercell%xcart,xred)
871 0 : xred_old = xred
872 0 : vel_cell(:,:) = zero
873 0 : vel(:,:) = zero
874 0 : gred(:,:) = zero
875 0 : fcart(:,:) = zero
876 :
877 : ! Run mover
878 : call mover(scfcv_args,ab_xfh,acell,effective_potential%crystal%amu,dtfil,electronpositron,&
879 : & rhog,rhor,dtset%rprimd_orig,vel,vel_cell,xred,xred_old,&
880 0 : & effective_potential=effective_potential,verbose=verbose,verbose_cycle=dtset%ionmov/=25,writeHIST=writeHIST)
881 :
882 0 : if(.not.effective_potential%anharmonics_terms%bounded)then
883 0 : write(message, '(2a)' ) ' => The model is not bounded'
884 : else
885 0 : write(message, '(2a)' ) ' => The model is bounded'
886 : end if
887 0 : call wrtout(std_out,message,'COLL')
888 : ! Exit if the model is bounded
889 0 : if(effective_potential%anharmonics_terms%bounded) then
890 0 : model_bound = jj
891 0 : model_ncoeffbound = ii
892 0 : exit
893 : end if
894 : end if
895 : end do
896 : end if
897 :
898 0 : ABI_FREE(list_tmp)
899 0 : ABI_FREE(list_bound)
900 0 : ABI_FREE(isPositive)
901 :
902 : ! Exit if the model is bounded
903 0 : if(effective_potential%anharmonics_terms%bounded) then
904 :
905 : ! Final transfert
906 0 : write(message, '(3a)' ) ch10,' => The model is now bounded'
907 0 : call wrtout(ab_out,message,'COLL')
908 0 : call wrtout(std_out,message,'COLL')
909 0 : do kk=ncoeff+1,ncoeff+model_ncoeffbound
910 0 : icoeff_bound = listcoeff_bound(model_bound,kk)-ncoeff ! need to remove ncoeff value
911 0 : call polynomial_coeff_free(coeffs_tmp(kk))
912 : call polynomial_coeff_init(coeff_values(model_bound,kk),&
913 : & coeffs_bound(icoeff_bound)%nterm,&
914 : & coeffs_tmp(kk),&
915 : & coeffs_bound(icoeff_bound)%terms,&
916 : & coeffs_bound(icoeff_bound)%name,&
917 0 : & check=.false.)
918 : end do
919 0 : ABI_FREE(coeff_values)
920 0 : ABI_FREE(listcoeff_bound)
921 0 : exit
922 : end if
923 0 : ABI_FREE(coeff_values)
924 :
925 : end do
926 :
927 0 : do ii=1,ncoeff_bound
928 0 : call polynomial_coeff_free(coeffs_bound(ii))
929 : end do
930 0 : ABI_SFREE(coeffs_bound)
931 :
932 : end if
933 :
934 0 : if(.not.effective_potential%anharmonics_terms%bounded)then
935 0 : write(message, '(3a)' ) ch10,' => The model cannot be bounded'
936 0 : call wrtout(ab_out,message,'COLL')
937 0 : call wrtout(std_out,message,'COLL')
938 0 : model_ncoeffbound = 0
939 0 : model_bound = 0
940 : end if
941 :
942 : ! Fit the final model
943 :
944 : call effective_potential_setCoeffs(coeffs_tmp(1:ncoeff+model_ncoeffbound),effective_potential,&
945 0 : & ncoeff+model_ncoeffbound)
946 :
947 : call fit_polynomial_coeff_fit(effective_potential,(/0/),(/0/),hist,0,(/0,0/),0,0,&
948 : & -1,inp%fit_nimposecoeff,inp%fit_imposecoeff,1,comm,verbose=.false.,fit_on=inp%fit_on,sel_on=inp%sel_on,&
949 0 : & max_nbody=inp%fit_max_nbody, drop_rate=0.0_dp, ncoeff_per_cycle=1, fit_weight_T=inp%fit_weight_T)
950 :
951 0 : write(message, '(3a)') ch10,' Fitted coefficients at the end of the fit bound process: '
952 0 : call wrtout(ab_out,message,'COLL')
953 0 : call wrtout(std_out,message,'COLL')
954 :
955 0 : do ii = 1,ncoeff+model_ncoeffbound
956 0 : write(message, '(a,I0,a,ES19.10,2a)') " ",ii," =>",&
957 0 : & effective_potential%anharmonics_terms%coefficients(ii)%coefficient,&
958 0 : & " ",trim(effective_potential%anharmonics_terms%coefficients(ii)%name)
959 0 : call wrtout(ab_out,message,'COLL')
960 0 : call wrtout(std_out,message,'COLL')
961 : end do
962 :
963 : ! Deallocation
964 0 : ABI_FREE(listcoeff)
965 0 : do ii=1,ncoeff_max
966 0 : call polynomial_coeff_free(coeffs_tmp(ii))
967 : end do
968 0 : ABI_SFREE(coeffs_tmp)
969 :
970 0 : do ii=1,ncoeff_max
971 0 : call polynomial_coeff_free(coeffs_all(ii))
972 : end do
973 0 : ABI_SFREE(coeffs_all)
974 :
975 : end if
976 :
977 : else if (option == -3) then
978 :
979 : !*************************************************************
980 : ! Call the routine for calculation of the energy for specific
981 : ! partern of displacement or strain for the effective
982 : ! Hamiltonian
983 : !*************************************************************
984 : ! write(message, '((80a),4a)' ) ('-',ii=1,80), ch10,&
985 : !& ' Effective Hamiltonian calculation'
986 : ! call wrtout(ab_out,message,'COLL')
987 : ! call wrtout(std_out,message,'COLL')
988 :
989 : ! acell = one
990 : ! call gstate(args_gs,acell,codvsn,cpui,dtfil,dtset,iexit,initialized,&
991 : !& mpi_enreg,npwtot,dtset%occ_orig,pawang,pawrad,pawtab,&
992 : !& psps,results_gs,dtset%rprimd_orig,scf_history,vel,vel_cell,wvl,xred)
993 :
994 : end if
995 :
996 : !***************************************************************
997 : ! 5 Deallocation of array
998 : !***************************************************************
999 :
1000 44 : ABI_FREE(gred)
1001 44 : ABI_FREE(fcart)
1002 44 : ABI_FREE(indsym)
1003 44 : ABI_FREE(rhog)
1004 44 : ABI_FREE(rhor)
1005 44 : ABI_FREE(vel)
1006 44 : ABI_FREE(xred)
1007 44 : ABI_FREE(xred_old)
1008 44 : ABI_FREE(ab_xfh%xfhist)
1009 :
1010 : ! if(option == -3)then
1011 : ! call args_gs_free(args_gs)
1012 : ! call psps_free(psps)
1013 : ! do ii = 1,npsp
1014 : ! call paw_setup_free(paw_setup(ii))
1015 : ! end do
1016 : ! ABI_FREE(paw_setup)
1017 : ! ABI_FREE(ipsp2xml)
1018 : ! ABI_FREE(pspheads)
1019 : ! call pawrad_free(pawrad)
1020 : ! call pawtab_free(pawtab)
1021 : ! ABI_FREE(pawrad)
1022 : ! ABI_FREE(pawtab)
1023 : ! ABI_FREE(npwtot)
1024 : ! end if
1025 44 : call dtset%free()
1026 44 : call results_gs%free()
1027 44 : call scfcv_destroy(scfcv_args)
1028 44 : call destroy_mpi_enreg(mpi_enreg)
1029 :
1030 : end if
1031 :
1032 3564 : write(message, '(a,(80a),a,a)' ) ch10,('=',ii=1,80),ch10
1033 44 : call wrtout(ab_out,message,'COLL')
1034 44 : call wrtout(std_out,message,'COLL')
1035 :
1036 88 : end subroutine mover_effpot
1037 : !!***
1038 :
1039 : !!****f* ABINIT/chksymmetrygroup
1040 : !! NAME
1041 : !! checksymmetrygroup
1042 : !!
1043 : !! FUNCTION
1044 : !! Find the ptgroup and symmetry relations (symre,tnons) of crystal by the lattice constants
1045 : !! rprimd and the reduced coordinates
1046 : !!
1047 : !! INPUTS
1048 : !! rprimd
1049 : !! xred
1050 : !! typat
1051 : !! msym : maximum symmetries defines sizes of symrel and tnons
1052 : !! natom
1053 : !!
1054 : !! OUTPUT
1055 : !! ptgroupma
1056 : !! spgroup: index of spacegroup
1057 : !! symrel(3,3,msym): symmetry relations
1058 : !! tnons(3,msym): translations
1059 : !!
1060 : !! SOURCE
1061 :
1062 1 : subroutine checksymmetrygroup(rprimd,xred,typat,msym,natom,ptgroupma,spgroup,symrel_out,tnons_out,nsym,tolsym)
1063 :
1064 : !Arguments ------------------------------------
1065 : !scalars
1066 : integer,intent(in) :: msym,natom
1067 : integer,intent(in) :: typat(natom)
1068 : integer,intent(out) :: ptgroupma,spgroup,nsym
1069 : real(dp),intent(inout) :: tolsym
1070 : ! Arrays
1071 : real(dp),intent(in) :: rprimd(3,3),xred(3,natom)
1072 : integer,intent(out) :: symrel_out(3,3,msym)
1073 : real(dp),intent(out) :: tnons_out(3,msym)
1074 :
1075 : !Local variables ---------------------------------------
1076 : !scalars
1077 : integer :: nptsym,use_inversion
1078 : integer :: chkprim
1079 : ! Arraiys
1080 2 : integer :: bravais(11),ptsymrel(3,3,msym)
1081 2 : integer :: symafm(msym),symrel(3,3,msym)
1082 2 : real(dp) :: gprimd(3,3),spinat(3,natom)
1083 1 : real(dp) :: tnons(3,msym)
1084 : real(dp) :: genafm(3)
1085 :
1086 : ! given the acel, rprim and coor
1087 : ! this suroutine find the symmetry group
1088 81 : spinat = 0
1089 1 : chkprim = 0
1090 1 : use_inversion = 0
1091 :
1092 : !write(std_out,*) "tolsym", tolsym, "tol3", tol3
1093 :
1094 1 : call symlatt(bravais,std_out,msym,nptsym,ptsymrel,rprimd,tol4)
1095 : !write(std_out,*) 'nptsym', nptsym
1096 1 : call matr3inv(rprimd,gprimd)
1097 : call symfind(gprimd,msym,natom,nptsym,0,nsym,&
1098 1 : & 0,ptsymrel,spinat,symafm,symrel,tnons,tolsym,typat,use_inversion,xred)
1099 :
1100 : !write(std_out,*) 'nsym', nsym
1101 1 : call symanal(bravais,chkprim,genafm,msym,nsym,ptgroupma,rprimd,spgroup,symafm,symrel,tnons,tol3)
1102 :
1103 : !write(std_out,*) 'nsym', nsym
1104 4993 : symrel_out = symrel
1105 1537 : tnons_out = tnons
1106 :
1107 44 : end subroutine checksymmetrygroup
1108 : !!***
1109 :
1110 : end module m_mover_effpot
1111 : !!***
|
