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