Branch data 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 : : !!***
|
