TABLE OF CONTENTS
ABINIT/erlxconv [ Functions ]
NAME
erlxconv
FUNCTION
FIXME: add description.
INPUTS
argin(sizein)=description
OUTPUT
argout(sizeout)=description
PARENTS
mover
CHILDREN
wrtout
SOURCE
1211 subroutine erlxconv(hist,iexit,itime,itime_hist,ntime,tolmxde) 1212 1213 !use m_abihist, only : abihist,abihist_findIndex 1214 1215 !This section has been created automatically by the script Abilint (TD). 1216 !Do not modify the following lines by hand. 1217 #undef ABI_FUNC 1218 #define ABI_FUNC 'erlxconv' 1219 use interfaces_14_hidewrite 1220 !End of the abilint section 1221 1222 implicit none 1223 1224 !Arguments ------------------------------------ 1225 !scalars 1226 integer,intent(in) :: itime,itime_hist,ntime 1227 integer,intent(inout) :: iexit 1228 real(dp), intent(in) :: tolmxde 1229 !arrays 1230 type(abihist),intent(inout) :: hist 1231 1232 !Local variables------------------------------- 1233 integer :: ihist,ihist_prev,ihist_prev2 1234 real(dp) :: ediff1,ediff2,maxediff 1235 character(len=500) :: message 1236 ! ************************************************************************* 1237 1238 if (itime_hist<3) then 1239 write(message, '(a,a,a)' ) ch10,& 1240 & ' erlxconv : minimum 3 Broyd/MD steps to check convergence of energy in relaxations',ch10 1241 call wrtout(std_out,message,'COLL') 1242 else 1243 ihist = hist%ihist 1244 ihist_prev = abihist_findIndex(hist,-1) 1245 ihist_prev2 = abihist_findIndex(hist,-2) 1246 ediff1 = hist%etot(ihist) - hist%etot(ihist_prev) 1247 ediff2 = hist%etot(ihist) - hist%etot(ihist_prev2) 1248 if ((abs(ediff1)<tolmxde).and.(abs(ediff2)<tolmxde)) then 1249 write(message, '(a,a,i4,a,a,a,a,a,es11.4,a,a)' ) ch10,& 1250 & ' At Broyd/MD step',itime,', energy is converged : ',ch10,& 1251 & ' the difference in energy with respect to the two ',ch10,& 1252 & ' previous steps is < tolmxde=',tolmxde,' ha',ch10 1253 call wrtout(ab_out,message,'COLL') 1254 call wrtout(std_out,message,'COLL') 1255 iexit=1 1256 else 1257 maxediff = max(abs(ediff1),abs(ediff2)) 1258 if(iexit==1)then 1259 write(message, '(a,a,a,a,i5,a,a,a,es11.4,a,es11.4,a,a)' ) ch10,& 1260 & ' erlxconv : WARNING -',ch10,& 1261 & ' ntime=',ntime,' was not enough Broyd/MD steps to converge energy: ',ch10,& 1262 & ' max difference in energy =',maxediff,' > tolmxde=',tolmxde,' ha',ch10 1263 call wrtout(std_out,message,'COLL') 1264 call wrtout(ab_out,message,'COLL') 1265 1266 write(std_out,"(8a)")ch10,& 1267 & "--- !RelaxConvergenceWarning",ch10,& 1268 & "message: | ",ch10,TRIM(indent(message)),ch10,& 1269 & "..." 1270 else 1271 write(message, '(a,a,i4,a,a,a,es11.4,a,es11.4,a,a)' ) ch10,& 1272 & ' erlxconv : at Broyd/MD step',itime,', energy has not converged yet. ',ch10,& 1273 & ' max difference in energy=',maxediff,' > tolmxde=',tolmxde,' ha',ch10 1274 call wrtout(std_out,message,'COLL') 1275 end if 1276 end if 1277 end if 1278 1279 end subroutine erlxconv
ABINIT/fconv [ Functions ]
NAME
fconv
FUNCTION
Check maximal absolute value of force (hartree/bohr) against
input tolerance; if below tolerance, return iexit=1.
Takes into account the fact that the Broyden (or moldyn) step
might be the last one (last itime), to print eventually modified message.
Stresses are also included in the check, provided that optcell/=0.
If optcell=1, takes only the trace into account
optcell=2, takes all components into account
optcell=3, takes traceless stress into account
optcell=4, takes sigma(1 1) into account
optcell=5, takes sigma(2 2) into account
optcell=6, takes sigma(3 3) into account
optcell=7, takes sigma(2,2),(2,3) and (3 3) into account
optcell=8, takes sigma(1,1),(1,3) and (3 3) into account
optcell=9, takes sigma(1,1),(1,2) and (2 2) into account
In the case of stresses, target the tensor strtarget, and
take into account the factor strfact
INPUTS
fcart(3,natom)= forces on atoms in hartree/bohr in cartesian coordinates iatfix(3,natom)=1 for frozen atom, 0 for unfrozen itime=current number of Broyden/Moldyn iterations natom=number of atoms in unit cell ntime=maximum number of Broyden/Moldyn iterations allowed optcell=option for taking stresses into account (see above) strfact=factor that multiplies the stresses when they are compared to forces. strtarget(6)=components of the target stress tensor (hartree/bohr^3) strten(6)=components of the stress tensor (hartree/bohr^3) tolmxf=tolerance on maximal absolute value of components of forces
OUTPUT
writes to unit std_out and to ab_out, and returns
SIDE EFFECTS
Input/Output at input : iexit= 0 if not the last itime, 1 if the last itime at output : iexit= 0 if not below tolerance, 1 if below tolerance
PARENTS
mover
CHILDREN
wrtout
SOURCE
1081 subroutine fconv(fcart,iatfix,iexit,itime,natom,ntime,optcell,strfact,strtarget,strten,tolmxf) 1082 1083 1084 !This section has been created automatically by the script Abilint (TD). 1085 !Do not modify the following lines by hand. 1086 #undef ABI_FUNC 1087 #define ABI_FUNC 'fconv' 1088 use interfaces_14_hidewrite 1089 !End of the abilint section 1090 1091 implicit none 1092 1093 !Arguments ------------------------------------ 1094 !scalars 1095 integer,intent(in) :: itime,natom,ntime,optcell 1096 integer,intent(inout) :: iexit 1097 real(dp),intent(in) :: strfact,tolmxf 1098 !arrays 1099 integer,intent(in) :: iatfix(3,natom) 1100 real(dp),intent(in) :: fcart(3,natom),strtarget(6),strten(6) 1101 1102 !Local variables------------------------------- 1103 !scalars 1104 integer :: iatom,idir,istr 1105 real(dp) :: fmax,strdiag 1106 character(len=500) :: message 1107 !arrays 1108 real(dp) :: dstr(6) 1109 1110 ! ************************************************************************* 1111 1112 !Compute maximal component of forces, EXCLUDING any fixed components 1113 fmax=zero 1114 do iatom=1,natom 1115 do idir=1,3 1116 if (iatfix(idir,iatom) /= 1) then 1117 if( abs(fcart(idir,iatom)) >= fmax ) fmax=abs(fcart(idir,iatom)) 1118 end if 1119 end do 1120 end do 1121 1122 dstr(:)=strten(:)-strtarget(:) 1123 1124 !Eventually take into account the stress 1125 if(optcell==1)then 1126 strdiag=(dstr(1)+dstr(2)+dstr(3))/3.0_dp 1127 if(abs(strdiag)*strfact >= fmax ) fmax=abs(strdiag)*strfact 1128 else if(optcell==2)then 1129 do istr=1,6 1130 if(abs(dstr(istr))*strfact >= fmax ) fmax=abs(dstr(istr))*strfact 1131 end do 1132 else if(optcell==3)then 1133 ! Must take away the trace from diagonal elements 1134 strdiag=(dstr(1)+dstr(2)+dstr(3))/3.0_dp 1135 do istr=1,3 1136 if(abs(dstr(istr)-strdiag)*strfact >= fmax ) fmax=abs(dstr(istr)-strdiag)*strfact 1137 end do 1138 do istr=4,6 1139 if(abs(dstr(istr))*strfact >= fmax ) fmax=abs(dstr(istr))*strfact 1140 end do 1141 else if(optcell==4 .or. optcell==5 .or. optcell==6)then 1142 if(abs(dstr(optcell-3))*strfact >= fmax ) fmax=abs(dstr(optcell-3))*strfact 1143 else if(optcell==7)then 1144 if(abs(dstr(2))*strfact >= fmax ) fmax=abs(dstr(2))*strfact 1145 if(abs(dstr(3))*strfact >= fmax ) fmax=abs(dstr(3))*strfact 1146 if(abs(dstr(4))*strfact >= fmax ) fmax=abs(dstr(4))*strfact 1147 else if(optcell==8)then 1148 if(abs(dstr(1))*strfact >= fmax ) fmax=abs(dstr(1))*strfact 1149 if(abs(dstr(3))*strfact >= fmax ) fmax=abs(dstr(3))*strfact 1150 if(abs(dstr(5))*strfact >= fmax ) fmax=abs(dstr(5))*strfact 1151 else if(optcell==9)then 1152 if(abs(dstr(1))*strfact >= fmax ) fmax=abs(dstr(1))*strfact 1153 if(abs(dstr(2))*strfact >= fmax ) fmax=abs(dstr(2))*strfact 1154 if(abs(dstr(6))*strfact >= fmax ) fmax=abs(dstr(6))*strfact 1155 end if 1156 1157 if (fmax<tolmxf) then 1158 write(message, '(a,a,i4,a,a,a,es11.4,a,es11.4,a,a)' ) ch10,& 1159 & ' At Broyd/MD step',itime,', gradients are converged : ',ch10,& 1160 & ' max grad (force/stress) =',fmax,' < tolmxf=',tolmxf,' ha/bohr (free atoms)',ch10 1161 call wrtout(ab_out,message,'COLL') 1162 call wrtout(std_out,message,'COLL') 1163 iexit=1 1164 else 1165 if(iexit==1)then 1166 write(message, '(a,a,a,a,i5,a,a,a,es11.4,a,es11.4,a,a)' ) ch10,& 1167 & ' fconv : WARNING -',ch10,& 1168 & ' ntime=',ntime,' was not enough Broyd/MD steps to converge gradients: ',ch10,& 1169 & ' max grad (force/stress) =',fmax,' > tolmxf=',tolmxf,' ha/bohr (free atoms)',ch10 1170 call wrtout(std_out,message,'COLL') 1171 call wrtout(ab_out,message,'COLL') 1172 1173 write(std_out,"(8a)")ch10,& 1174 & "--- !RelaxConvergenceWarning",ch10,& 1175 & "message: | ",ch10,TRIM(indent(message)),ch10,& 1176 & "..." 1177 1178 else 1179 write(message, '(a,i4,a,a,a,es11.4,a,es11.4,a,a)' ) & 1180 & ' fconv : at Broyd/MD step',itime,', gradients have not converged yet. ',ch10,& 1181 & ' max grad (force/stress) =',fmax,' > tolmxf=',tolmxf,' ha/bohr (free atoms)',ch10 1182 call wrtout(std_out,message,'COLL') 1183 end if 1184 iexit=0 1185 end if 1186 1187 end subroutine fconv
ABINIT/mover [ Functions ]
NAME
mover
FUNCTION
Move ion or change acell acording to forces and stresses
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, SE) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
amass(natom)=mass of each atom, in unit of electronic mass (=amu*1822...) dtfil <type(datafiles_type)>=variables related to files dtset <type(dataset_type)>=all input variables for this dataset | mband=maximum number of bands | mgfft=maximum size of 1D FFTs | mkmem =number of k points treated by this node | angular momentum for nonlocal pseudopotential | mpw=maximum dimensioned size of npw. | natom=number of atoms in unit cell | except on first call (hartree/bohr); updated on output | nfft=(effective) number of FFT grid points (for this processor) | for the "coarse" grid (see NOTES below) | nkpt=number of k points. | nspden=number of spin-density components | nsppol=1 for unpolarized, 2 for spin-polarized | nsym=number of symmetry elements in space group mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol mpi_enreg=informations about MPI parallelization nfftf=(effective) number of FFT grid points (for this processor) for the "fine" grid (see NOTES below) npwarr(nkpt)=number of planewaves in basis and boundary at this k point. nattyp(ntypat)= # atoms of each type. paw_dmft <type(paw_dmft_type)>= paw+dmft related data psps <type(pseudopotential_type)>=variables related to pseudopotentials | mpsang= 1+maximum angular momentum for nonlocal pseudopotentials rprimd(3,3)=dimensional primitive translations (bohr)
OUTPUT
results_gs <type(results_gs_type)>=results (energy and its components, forces and its components, the stress tensor) of a ground-state computation eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree) resid(mband*nkpt*nsppol)=residuals for each band over all k points.
SIDE EFFECTS
Rest of i/o is related to lda acell(3)=length scales of primitive translations (bohr) cg(2,mcg)=array for planewave coefficients of wavefunctions. electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation initialized= if 0 the initialisation of the gstate run is not yet finished irrzon(nfft**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data pawrhoij(natom*usepaw) <type(pawrhoij_type)>= -PAW only- atomic occupancies occ(mband*nkpt*nsppol=occupation number for each band (usually 2) at each k point. rhog(2,nfftf)=array for Fourier transform of electron density rhor(nfftf,nspden)=array for electron density in electrons/bohr**3. scf_history <type(scf_history_type)>=arrays obtained from previous SCF cycles symrec(3,3,nsym)=symmetry operations in reciprocal space taug(2,nfftf*dtset%usekden)=array for Fourier transform of kinetic energy density taur(nfftf,nspden*dtset%usekden)=array for kinetic energy density vel(3,natom)=old value of velocity; updated on output vel_cell(3,3)=old value of cell parameters velocity; updated on output xred(3,natom)=reduced dimensionless atomic coordinates; updated on output xred_old(3,natom)=work space for old xred eff_pot<type(effective_potential_type)> = optional,effective_potential datatype verbose = optional, default is true, flag to disable the verbose mode write_HIST = optional, default is true, flag to disble the write of the HIST file
NOTES
This subroutine uses the arguments natom, xred, vel, amass, vis, and dtion (the last two contained in dtset) to make molecular dynamics updates. The rest of the lengthy argument list supports the underlying lda computation of forces, returned from subroutine scfcv USE OF FFT GRIDS: ================= In case of PAW: --------------- Two FFT grids are used: - A "coarse" FFT grid (defined by ecut) for the application of the Hamiltonian on the plane waves basis. It is defined by nfft, ngfft, mgfft, ... Hamiltonian, wave-functions, density related to WFs (rhor here), ... are expressed on this grid. - A "fine" FFT grid (defined) by ecutdg) for the computation of the density inside PAW spheres. It is defined by nfftf, ngfftf, mgfftf, ... Total density, potentials, ... are expressed on this grid. In case of norm-conserving: --------------------------- - Only the usual FFT grid (defined by ecut) is used. It is defined by nfft, ngfft, mgfft, ... For compatibility reasons, (nfftf,ngfftf,mgfftf) are set equal to (nfft,ngfft,mgfft) in that case.
PARENTS
gstate,mover_effpot
CHILDREN
abiforstr_fin,abiforstr_ini,abihist_bcast,abihist_compare_and_copy abihist_free,abihist_init,abimover_fin,abimover_ini,abimover_nullify chkdilatmx,crystal_free,crystal_init,dtfil_init_time effective_potential_evaluate,erlxconv,fcart2fred,fconv,hist2var initylmg,matr3inv,mttk_fin,mttk_ini,prec_simple,pred_bfgs,pred_delocint pred_diisrelax,pred_hmc,pred_isokinetic,pred_isothermal,pred_langevin pred_lbfgs,pred_lotf,pred_moldyn,pred_nose,pred_simple,pred_srkna14 pred_steepdesc,pred_velverlet,pred_verlet,prtxfase,read_md_hist scfcv_run,status,symmetrize_xred,var2hist,vel2hist,write_md_hist wrt_moldyn_netcdf,wrtout,wvl_mkrho,wvl_wfsinp_reformat,xfh_update xmpi_barrier,xmpi_isum,xmpi_wait
SOURCE
119 #if defined HAVE_CONFIG_H 120 #include "config.h" 121 #endif 122 123 #include "abi_common.h" 124 125 subroutine mover(scfcv_args,ab_xfh,acell,amass,dtfil,& 126 & electronpositron,rhog,rhor,rprimd,vel,vel_cell,xred,xred_old,& 127 & effective_potential,verbose,writeHIST) 128 129 use defs_basis 130 use defs_abitypes 131 use m_errors 132 use m_abimover 133 use m_abihist 134 use m_xmpi 135 use m_nctk 136 #ifdef HAVE_NETCDF 137 use netcdf 138 #endif 139 #if defined HAVE_LOTF 140 use lotfpath 141 use m_pred_lotf 142 #endif 143 144 use m_fstrings, only : strcat, sjoin, indent 145 use m_geometry, only : fcart2fred, chkdilatmx 146 use m_crystal, only : crystal_init, crystal_free, crystal_t 147 use m_crystal_io, only : crystal_ncwrite_path 148 use m_time, only : abi_wtime, sec2str 149 use m_exit, only : get_start_time, have_timelimit_in, get_timelimit, enable_timelimit_in 150 use m_electronpositron, only : electronpositron_type 151 use m_scfcv, only : scfcv_t, scfcv_run 152 use m_effective_potential,only : effective_potential_type,effective_potential_evaluate 153 154 !This section has been created automatically by the script Abilint (TD). 155 !Do not modify the following lines by hand. 156 #undef ABI_FUNC 157 #define ABI_FUNC 'mover' 158 use interfaces_14_hidewrite 159 use interfaces_32_util 160 use interfaces_41_geometry 161 use interfaces_45_geomoptim 162 use interfaces_56_recipspace 163 use interfaces_59_ionetcdf 164 use interfaces_67_common 165 use interfaces_79_seqpar_mpi 166 use interfaces_95_drive, except_this_one => mover 167 !End of the abilint section 168 169 implicit none 170 171 !Arguments ------------------------------------ 172 !scalars 173 type(scfcv_t),intent(inout) :: scfcv_args 174 type(datafiles_type),intent(inout),target :: dtfil 175 type(electronpositron_type),pointer :: electronpositron 176 type(ab_xfh_type),intent(inout) :: ab_xfh 177 type(effective_potential_type),optional,intent(inout) :: effective_potential 178 logical,optional,intent(in) :: verbose 179 logical,optional,intent(in) :: writeHIST 180 !arrays 181 real(dp),intent(inout) :: acell(3) 182 real(dp), intent(in),target :: amass(:) !(scfcv%dtset%natom) cause segfault of g95 on yquem_g95 A check of the dim has been added 183 real(dp), pointer :: rhog(:,:),rhor(:,:) 184 real(dp), intent(inout) :: xred(3,scfcv_args%dtset%natom),xred_old(3,scfcv_args%dtset%natom) 185 real(dp), intent(inout) :: vel(3,scfcv_args%dtset%natom),vel_cell(3,3),rprimd(3,3) 186 187 !Local variables------------------------------- 188 !scalars 189 integer,parameter :: level=102,master=0 190 type(abihist) :: hist,hist_prev 191 type(abimover) :: ab_mover 192 type(abimover_specs) :: specs 193 type(abiforstr) :: preconforstr ! Preconditioned forces and stress 194 type(mttk_type) :: mttk_vars 195 integer :: itime,icycle,itime_hist,iexit=0,ifirst,ihist_prev,ihist_prev2,timelimit_exit,ncycle,nhisttot,kk,jj,me 196 integer :: nloop,ntime,option,comm 197 integer :: nerr_dilatmx,my_quit,ierr,quitsum_request 198 integer ABI_ASYNC :: quitsum_async 199 character(len=500) :: message 200 character(len=500) :: dilatmx_errmsg 201 character(len=8) :: stat4xml 202 character(len=35) :: fmt 203 character(len=fnlen) :: filename 204 real(dp) :: favg 205 logical :: DEBUG=.FALSE., need_verbose=.TRUE.,need_writeHIST=.TRUE. 206 logical :: need_scfcv_cycle = .TRUE. 207 logical :: change,useprtxfase 208 logical :: skipcycle 209 integer :: minIndex,ii,similar,conv_retcode 210 integer :: iapp 211 real(dp) :: minE,wtime_step,now,prev 212 type(crystal_t) :: crystal 213 !arrays 214 real(dp) :: gprimd(3,3),rprim(3,3),rprimd_prev(3,3) 215 real(dp),allocatable :: amu(:),fred_corrected(:,:),xred_prev(:,:) 216 217 ! *************************************************************** 218 need_verbose=.TRUE. 219 if(present(verbose)) need_verbose = verbose 220 221 need_writeHIST=.TRUE. 222 if(present(writeHIST)) need_writeHIST = writeHIST 223 224 call status(0,dtfil%filstat,iexit,level,'init ') 225 226 if ( size(amass) /= scfcv_args%dtset%natom ) then 227 MSG_BUG("amass does not have the proper size.") 228 end if 229 230 ! enable time limit handler if not done in callers. 231 if (need_verbose .and. enable_timelimit_in(ABI_FUNC) == ABI_FUNC) then 232 write(std_out,*)"Enabling timelimit check in function: ",trim(ABI_FUNC)," with timelimit: ",trim(sec2str(get_timelimit())) 233 end if 234 235 !Table of contents 236 !(=>) Refers to an important call (scfcv,pred_*) 237 ! 238 !01. Initialization of indexes and allocations of arrays 239 !02. Particularities of each predictor 240 !03. Set the number of iterations ntime 241 !04. Try to read history of previous calculations 242 !05. Allocate the hist structure 243 !06. First output before any itime or icycle 244 !07. Fill the history of the first SCFCV 245 !08. Loop for itime (From 1 to ntime) 246 !09. Loop for icycle (From 1 to ncycles) 247 !10. Output for each icycle (and itime) 248 !11. Symmetrize atomic coordinates over space group elements 249 !12. => Call to SCFCV routine and fill history with forces 250 !13. Write the history into the _HIST file 251 !14. Output after SCFCV 252 !15. => Test Convergence of forces and stresses 253 !16. => Precondition forces, stress and energy 254 !17. => Call to each predictor 255 !18. Use the history to extract the new values 256 !19. End loop icycle 257 !20. End loop itime 258 !21. Set the final values of xred 259 !22. XML Output at the end 260 !23. Deallocate hist and ab_mover datatypes 261 ! 262 !IONMOV values: 263 ! 264 !1. Molecular dynamics without viscosity (vis=0) 265 !1. Molecular dynamics with viscosity (vis/=0) 266 !2. Broyden-Fletcher-Goldfard-Shanno method (forces) 267 !3. Broyden-Fletcher-Goldfard-Shanno method (forces,Tot energy) 268 !4. Conjugate gradient of potential and ionic degrees of freedom 269 !5. Simple relaxation of ionic positions 270 !6. Verlet algorithm for molecular dynamics 271 !7. Verlet algorithm blocking every atom where dot(vel,force)<0 272 !8. Verlet algorithm with a nose-hoover thermostat 273 !9. Langevin molecular dynamics 274 !10. BFGS with delocalized internal coordinates 275 !11. Conjugate gradient algorithm 276 !12. Isokinetic ensemble molecular dynamics 277 !13. Isothermal/isenthalpic ensemble molecular dynamics 278 !14. Symplectic algorithm Runge-Kutta-Nystrom SRKNa14 279 !20. Ionic positions relaxation using DIIS 280 !21. Steepest descent algorithm 281 !23. LOTF method 282 !24. Velocity verlet molecular dynamics 283 284 call abimover_nullify(ab_mover) 285 286 call abimover_ini(ab_mover,specs,& 287 & scfcv_args%dtset%delayperm,& 288 & scfcv_args%dtset%diismemory,& 289 & scfcv_args%dtset%goprecon,& 290 & scfcv_args%dtset%jellslab,& 291 & scfcv_args%dtset%natom,& 292 & scfcv_args%dtset%nconeq,& 293 & scfcv_args%dtset%nnos,& 294 & scfcv_args%dtset%nsym,& 295 & scfcv_args%dtset%ntypat,& 296 & scfcv_args%dtset%optcell,& 297 & scfcv_args%dtset%restartxf,& 298 & scfcv_args%dtset%signperm,& 299 & scfcv_args%dtset%ionmov,& 300 & scfcv_args%dtset%bmass,& 301 & scfcv_args%dtset%dtion,& 302 & scfcv_args%dtset%friction,& 303 & scfcv_args%dtset%mdwall,& 304 & scfcv_args%dtset%noseinert,& 305 & scfcv_args%dtset%strprecon,& 306 & scfcv_args%dtset%vis,& 307 & scfcv_args%dtset%iatfix,& 308 & scfcv_args%dtset%symrel,& 309 & scfcv_args%dtset%typat,& 310 & scfcv_args%dtset%prtatlist,& 311 & amass,& 312 & scfcv_args%dtset%goprecprm,& 313 & scfcv_args%dtset%mdtemp,& 314 & scfcv_args%dtset%strtarget,& 315 & scfcv_args%dtset%qmass,& 316 & scfcv_args%dtset%znucl,& 317 & dtfil%fnameabi_hes,& 318 & dtfil%filnam_ds) 319 320 if (ab_mover%ionmov==13)then 321 call mttk_ini(mttk_vars,ab_mover%nnos) 322 end if 323 324 !########################################################### 325 !### 03. Set the number of iterations ntime 326 !### By default ntime==1 but if the user enter a lower 327 !### value mover will execute at least one iteration 328 329 if (scfcv_args%dtset%ntime<1)then 330 ntime=1 331 else 332 ntime=scfcv_args%dtset%ntime 333 end if 334 335 !########################################################### 336 !### 04. Try to read history of previous calculations 337 !### It requires access to the NetCDF library 338 339 !Init MPI data 340 comm=scfcv_args%mpi_enreg%comm_cell 341 me=xmpi_comm_rank(comm) 342 343 #if defined HAVE_NETCDF 344 filename=trim(ab_mover%filnam_ds(4))//'_HIST.nc' 345 346 if (ab_mover%restartxf<0)then 347 ! Read history from file (and broadcast if MPI) 348 if (me==master) then 349 call read_md_hist(filename,hist_prev,specs%isVused,specs%isARused,ab_mover%restartxf==-3) 350 end if 351 call abihist_bcast(hist_prev,master,comm) 352 353 ! If restartxf specifies to reconstruct the history 354 if (hist_prev%mxhist>0.and.ab_mover%restartxf==-1)then 355 ntime=ntime+hist_prev%mxhist 356 end if 357 358 ! If restartxf specifies to start from the lowest energy 359 if (hist_prev%mxhist>0.and.ab_mover%restartxf==-2)then 360 minE=hist_prev%etot(1) 361 minIndex=1 362 do ii=1,hist_prev%mxhist 363 if(need_verbose) write(std_out,*) 'Iteration:',ii,' Total Energy:',hist_prev%etot(ii) 364 if (minE>hist_prev%etot(ii))then 365 minE=hist_prev%etot(ii) 366 minIndex=ii 367 end if 368 end do 369 if(need_verbose)write(std_out,*) 'The lowest energy occurs at iteration:',minIndex,'etotal=',minE 370 acell(:) =hist_prev%acell(:,minIndex) 371 rprimd(:,:)=hist_prev%rprimd(:,:,minIndex) 372 xred(:,:) =hist_prev%xred(:,:,minIndex) 373 call abihist_free(hist_prev) 374 end if 375 ! If restarxf specifies to start to the last iteration 376 if (hist_prev%mxhist>0.and.ab_mover%restartxf==-3)then 377 acell(:) =hist_prev%acell(:,hist_prev%mxhist) 378 rprimd(:,:)=hist_prev%rprimd(:,:,hist_prev%mxhist) 379 xred(:,:) =hist_prev%xred(:,:,hist_prev%mxhist) 380 call abihist_free(hist_prev) 381 end if 382 383 end if !if (ab_mover%restartxf<=0) 384 385 #endif 386 387 !########################################################### 388 !### 05. Allocate the hist structure 389 390 iexit=0; timelimit_exit=0 391 ncycle=specs%ncycle 392 393 nhisttot=ncycle*ntime;if (scfcv_args%dtset%nctime>0) nhisttot=nhisttot+1 394 395 !AM_2017 New version of the hist, we just store the needed history step not all of them... 396 if(specs%nhist/=-1) nhisttot = specs%nhist ! We don't need to store all the history 397 call abihist_init(hist,ab_mover%natom,nhisttot,specs%isVused,specs%isARused) 398 call abiforstr_ini(preconforstr,ab_mover%natom) 399 400 !########################################################### 401 !### 06. First output before any itime or icycle 402 403 !If effective potential is present, 404 ! forces will be compute with it 405 if (present(effective_potential)) then 406 need_scfcv_cycle = .FALSE. 407 if(need_verbose)then 408 write(message,'(2a,i2,5a,80a)')& 409 & ch10,'=== [ionmov=',ab_mover%ionmov,'] ',trim(specs%method),' with effective potential',& 410 & ch10,('=',kk=1,80) 411 call wrtout(ab_out,message,'COLL') 412 call wrtout(std_out,message,'COLL') 413 end if 414 else 415 if(need_verbose)then 416 write(message,'(a,a,i2,a,a,a,80a)')& 417 & ch10,'=== [ionmov=',ab_mover%ionmov,'] ',specs%method,& 418 & ch10,('=',kk=1,80) 419 call wrtout(ab_out,message,'COLL') 420 call wrtout(std_out,message,'COLL') 421 end if 422 end if 423 424 !Format for printing on each cycle 425 write(fmt,'(a6,i2,a4,i2,a4,i2,a4,i2,a9)')& 426 & '(a,a,i',int(log10(real(ntime))+1),& 427 & ',a,i',int(log10(real(ntime))+1),& 428 & ',a,i',int(log10(real(ncycle))+1),& 429 & ',a,i',int(log10(real(ncycle))+1),& 430 & ',a,a,80a)' 431 432 !########################################################### 433 !### 07. Fill the history of the first SCFCV 434 435 !Compute atomic masses (in a.u.) 436 ABI_ALLOCATE(amu,(scfcv_args%dtset%ntypat)) 437 do kk=1,ab_mover%ntypat 438 do jj=1,ab_mover%natom 439 if (kk==ab_mover%typat(jj)) then 440 amu(kk)=amass(jj)/amu_emass 441 exit 442 end if 443 end do 444 end do 445 446 !Fill history with the values of xred, acell and rprimd 447 call var2hist(acell,hist,ab_mover%natom,rprimd,xred,DEBUG) 448 449 !Fill velocities and ionic kinetic energy 450 call vel2hist(ab_mover%amass,hist,vel,vel_cell) 451 hist%time(hist%ihist)=zero 452 453 !Decide if prtxfase will be called 454 useprtxfase=.FALSE. 455 do ii=1,ab_mover%natom 456 if (ab_mover%prtatlist(ii)/=0)then 457 useprtxfase=.TRUE. 458 exit 459 end if 460 end do 461 462 !At beginning no error 463 nerr_dilatmx = 0 464 465 ABI_ALLOCATE(xred_prev,(3,scfcv_args%dtset%natom)) 466 467 !########################################################### 468 !### 08. Loop for itime (From 1 to ntime) 469 quitsum_request = xmpi_request_null 470 471 do itime=1,ntime 472 473 ! Handle time limit condition. 474 if (itime == 1) prev = abi_wtime() 475 if (itime > 1) then 476 now = abi_wtime() 477 wtime_step = now - prev 478 prev = now 479 write(message,*)sjoin("mover: previous time step took ",sec2str(wtime_step)) 480 if(need_verbose)call wrtout(std_out, message, "COLL") 481 if (have_timelimit_in(ABI_FUNC)) then 482 if (itime > 2) then 483 call xmpi_wait(quitsum_request,ierr) 484 if (quitsum_async > 0) then 485 write(message,"(3a)")"Approaching time limit ",trim(sec2str(get_timelimit())),& 486 & ". Will exit itime loop in mover." 487 if(need_verbose)MSG_COMMENT(message) 488 if(need_verbose)call wrtout(ab_out, message, "COLL") 489 timelimit_exit = 1 490 exit 491 end if 492 end if 493 494 my_quit = 0; if (now - get_start_time() + 2.15 * wtime_step > get_timelimit()) my_quit = 1 495 call xmpi_isum(my_quit,quitsum_async,comm,quitsum_request,ierr) 496 end if 497 end if 498 499 skipcycle=.FALSE. 500 #if defined HAVE_LOTF 501 if(ab_mover%ionmov==23 .and. .not. lotf_extrapolation(itime)) skipcycle=.True. 502 #endif 503 504 ! ########################################################### 505 ! ### 09. Loop for icycle (From 1 to ncycles) 506 do icycle=1,ncycle 507 itime_hist = (itime-1)*ncycle + icycle ! Store the time step in of the history 508 509 ! ########################################################### 510 ! ### 10. Output for each icycle (and itime) 511 if(need_verbose)then 512 write(message,fmt)& 513 & ch10,'--- Iteration: (',itime,'/',ntime,') Internal Cycle: (',icycle,'/',ncycle,')',ch10,('-',kk=1,80) 514 call wrtout(ab_out,message,'COLL') 515 call wrtout(std_out,message,'COLL') 516 end if 517 if (useprtxfase) then 518 call prtxfase(ab_mover,hist,itime_hist,std_out,mover_BEFORE) 519 end if 520 521 xred_prev(:,:)=xred(:,:) 522 rprimd_prev(:,:)=rprimd(:,:) 523 524 ! ########################################################### 525 ! ### 11. Symmetrize atomic coordinates over space group elements 526 527 call symmetrize_xred(scfcv_args%indsym,ab_mover%natom,& 528 & scfcv_args%dtset%nsym,scfcv_args%dtset%symrel,scfcv_args%dtset%tnons,xred) 529 530 change=any(xred(:,:)/=xred_prev(:,:)) 531 if (change)then 532 hist%xred(:,:,hist%ihist)=xred(:,:) 533 if(need_verbose) then 534 write(std_out,*) 'WARNING: ATOMIC COORDINATES WERE SYMMETRIZED' 535 write(std_out,*) 'DIFFERENCES:' 536 do kk=1,ab_mover%natom 537 write(std_out,*) xred(:,kk)-xred_prev(:,kk) 538 end do 539 end if 540 xred_prev(:,:)=xred(:,:) 541 end if 542 543 ! ########################################################### 544 ! ### 12. => Call to SCFCV routine and fill history with forces 545 if (need_verbose) then 546 if (need_scfcv_cycle) then 547 write(message,'(a,3a,33a,44a)')& 548 & ch10,('-',kk=1,3),& 549 & 'SELF-CONSISTENT-FIELD CONVERGENCE',('-',kk=1,44) 550 else 551 write(message,'(a,3a,33a,44a)')& 552 & ch10,('-',kk=1,3),& 553 & 'EFFECTIVE POTENTIAL CALCULATION',('-',kk=1,44) 554 end if 555 call wrtout(ab_out,message,'COLL') 556 call wrtout(std_out,message,'COLL') 557 end if 558 559 if(hist_prev%mxhist>0.and.ab_mover%restartxf==-1.and.hist_prev%ihist<=hist_prev%mxhist)then 560 561 call abihist_compare_and_copy(hist_prev,hist,ab_mover%natom,similar,tol8,specs%nhist==nhisttot) 562 hist_prev%ihist=hist_prev%ihist+1 563 564 else 565 scfcv_args%ndtpawuj=0 566 iapp=itime 567 if(icycle>1.and.icycle/=ncycle) iapp=-1 568 if(itime==1 .and. icycle/=ncycle ) iapp=-icycle-1 569 if (ab_mover%ionmov==14.and.(icycle<ncycle)) iapp=-1 570 571 #if defined HAVE_LOTF 572 if (ab_mover%ionmov/=23 .or.(lotf_extrapolation(itime).and.(icycle/=1.or.itime==1)))then 573 #endif 574 !call scfcv_new2(scfcv_args,electronpositron,rhog,rhor,rprimd,xred,xred_old,conv_retcode) 575 576 !WVL - reformat the wavefunctions in the case of xred != xred_old 577 if (scfcv_args%dtset%usewvl == 1 .and. maxval(xred_old - xred) > zero) then 578 !Before running scfcv, on non-first geometry step iterations, 579 ! we need to reformat the wavefunctions, taking into acount the new 580 ! coordinates. We prepare to change rhog (to be removed) and rhor. 581 ABI_DEALLOCATE(rhog) 582 ABI_DEALLOCATE(rhor) 583 call wvl_wfsinp_reformat(scfcv_args%dtset, scfcv_args%mpi_enreg,& 584 & scfcv_args%psps, rprimd, scfcv_args%wvl, xred, xred_old) 585 scfcv_args%nfftf = scfcv_args%dtset%nfft 586 ABI_ALLOCATE(rhog,(2, scfcv_args%dtset%nfft)) 587 ABI_ALLOCATE(rhor,(2, scfcv_args%dtset%nfft)) 588 call wvl_mkrho(scfcv_args%dtset, scfcv_args%irrzon, scfcv_args%mpi_enreg,& 589 & scfcv_args%phnons, rhor,scfcv_args%wvl%wfs,scfcv_args%wvl%den) 590 end if 591 592 ! MAIN CALL TO SELF-CONSISTENT FIELD ROUTINE 593 if (need_scfcv_cycle) then 594 call dtfil_init_time(dtfil,iapp) 595 call scfcv_run(scfcv_args,electronpositron,rhog,rhor,rprimd,xred,xred_old,conv_retcode) 596 if (conv_retcode == -1) then 597 message = "Scf cycle returned conv_retcode == -1 (timelimit is approaching), this should not happen inside mover" 598 MSG_WARNING(message) 599 end if 600 601 else 602 ! For monte carlo don't need to recompute energy here 603 ! (done in pred_montecarlo) 604 call effective_potential_evaluate( & 605 & effective_potential,scfcv_args%results_gs%etotal,& 606 & scfcv_args%results_gs%fcart,scfcv_args%results_gs%fred,& 607 & scfcv_args%results_gs%strten,ab_mover%natom,rprimd,xred=xred,verbose=need_verbose) 608 609 ! Check if the simulation does not diverged... 610 if(itime > 3 .and.ABS(scfcv_args%results_gs%etotal - hist%etot(1)) > 1E4)then 611 ! We set to false the flag corresponding to the bound 612 effective_potential%anharmonics_terms%bounded = .FALSE. 613 if(need_verbose.and.me==master)then 614 message = "The simulation is diverging, please check your effective potential" 615 MSG_WARNING(message) 616 end if 617 ! Set the flag to finish the simulation 618 iexit=1 619 stat4xml="Failed" 620 else 621 ! We set to true the flag corresponding to the bound 622 effective_potential%anharmonics_terms%bounded = .TRUE. 623 end if 624 end if 625 #if defined HAVE_LOTF 626 end if 627 #endif 628 ! ANOMALOUS SITUATION 629 ! This is the only case where rprimd could change inside scfcv 630 ! It generates an weird condition, we start with a certain 631 ! value for rprimd before scfcv and after we finish with 632 ! a different value. 633 ! Notice that normally scfcv should not change rprimd 634 ! And even worse if optcell==0 635 ! The solution here is to recompute acell and store these value 636 ! in the present record even if initially it was not exactly 637 ! the value entering in scfcv 638 ! One test case with these condition is bigdft/t10 639 if (any(rprimd(:,:)/=rprimd_prev(:,:))) then 640 hist%acell(:,hist%ihist)=acell(:) 641 hist%rprimd(:,:,hist%ihist)=rprimd(:,:) 642 end if 643 644 ! ANOMALOUS SITUATIONS 645 ! * In ionmov 4 & 5 xred could change inside SCFCV 646 ! So we need to take the values from the output 647 ! 648 ! * Inside scfcv.F90 there is a call to symmetrize_xred.F90 649 ! for the first SCF cycle symmetrize_xred could change xred 650 if (ab_mover%ionmov<10)then 651 change=any(xred(:,:)/=xred_prev(:,:)) 652 if (change)then 653 hist%xred(:,:,hist%ihist)=xred(:,:) 654 if(need_verbose)then 655 write(std_out,*) 'WARNING: ATOMIC COORDINATES WERE SYMMETRIZED AFTER SCFCV' 656 write(std_out,*) 'DIFFERENCES:' 657 do kk=1,ab_mover%natom 658 write(std_out,*) xred(:,kk)-xred_prev(:,kk) 659 end do 660 end if 661 end if !if (change) 662 end if 663 664 ! Fill velocities and ionic kinetic energy 665 call vel2hist(ab_mover%amass,hist,vel,vel_cell) 666 667 hist%fcart(:,:,hist%ihist)=scfcv_args%results_gs%fcart(:,:) 668 hist%strten(:,hist%ihist) =scfcv_args%results_gs%strten(:) 669 hist%etot(hist%ihist) =scfcv_args%results_gs%etotal 670 hist%entropy(hist%ihist) =scfcv_args%results_gs%energies%entropy 671 hist%time(hist%ihist) =real(itime,kind=dp) 672 673 ! !###################################################################### 674 end if ! if(hist_prev%mxhist>0.and.ab_mover%restartxf==-1.and.hist_prev%ihist<=hist_prev%mxhist)then 675 676 ! Store trajectory in xfh file 677 if((ab_xfh%nxfh==0.or.itime/=1)) then 678 ABI_ALLOCATE(fred_corrected,(3,scfcv_args%dtset%natom)) 679 call fcart2fred(hist%fcart(:,:,hist%ihist),fred_corrected,rprimd,ab_mover%natom) 680 ! Get rid of mean force on whole unit cell, 681 ! but only if no generalized constraints are in effect 682 if (ab_mover%nconeq==0)then 683 do ii=1,3 684 if (ii/=3.or.ab_mover%jellslab==0) then 685 favg=sum(fred_corrected(ii,:))/dble(ab_mover%natom) 686 fred_corrected(ii,:)=fred_corrected(ii,:)-favg 687 end if 688 end do 689 end if 690 if (ncycle<10.and.ab_mover%restartxf>=0) then 691 do ii=1,3 692 rprim(ii,1:3)=rprimd(ii,1:3)/acell(1:3) 693 end do 694 695 ! The size of ab_xfh%xfhist is to big for very large supercell. 696 ! Call it only for specific ionmov 697 if(any((/2,3,10,11,22/)==ab_mover%ionmov)) then 698 call xfh_update(ab_xfh,acell,fred_corrected,ab_mover%natom,rprim,& 699 & hist%strten(:,hist%ihist),xred) 700 end if 701 end if 702 ABI_DEALLOCATE(fred_corrected) 703 end if 704 705 ! ########################################################### 706 ! ### 13. Write the history into the _HIST file 707 ! ### 708 709 #if defined HAVE_NETCDF 710 if (need_writeHIST.and.me==master) then 711 ifirst=merge(0,1,(itime>1.or.icycle>1)) 712 call write_md_hist(hist,filename,ifirst,itime_hist,ab_mover%natom,scfcv_args%dtset%nctime,& 713 & ab_mover%ntypat,ab_mover%typat,amu,ab_mover%znucl,ab_mover%dtion,scfcv_args%dtset%mdtemp) 714 end if 715 #endif 716 717 ! ########################################################### 718 ! ### 14. Output after SCFCV 719 if(need_verbose.and.need_scfcv_cycle)then 720 write(message,'(a,3a,a,72a)')& 721 & ch10,('-',kk=1,3),'OUTPUT',('-',kk=1,71) 722 call wrtout(ab_out,message,'COLL') 723 call wrtout(std_out,message,'COLL') 724 end if 725 if (useprtxfase) then 726 call prtxfase(ab_mover,hist,itime_hist,ab_out,mover_AFTER) 727 call prtxfase(ab_mover,hist,itime_hist,std_out,mover_AFTER) 728 end if 729 730 ! ########################################################### 731 ! ### 15. => Test Convergence of forces and stresses 732 733 if (itime==ntime.and.icycle==ncycle)then 734 iexit=1 735 stat4xml="Failed" 736 else 737 stat4xml="Succeded" 738 end if 739 740 ! Only if convergence is needed 741 if(specs%isFconv)then 742 if ((ab_mover%ionmov/=4.and.ab_mover%ionmov/=5).or.mod(itime,2)==1)then 743 if (scfcv_args%dtset%tolmxf/=0)then 744 call fconv(hist%fcart(:,:,hist%ihist),& 745 & scfcv_args%dtset%iatfix, & 746 & iexit,itime,& 747 & ab_mover%natom,& 748 & ntime,& 749 & ab_mover%optcell,& 750 & scfcv_args%dtset%strfact,& 751 & scfcv_args%dtset%strtarget,& 752 & hist%strten(:,hist%ihist),& 753 & scfcv_args%dtset%tolmxf) 754 else 755 call erlxconv(hist,iexit,itime,itime_hist,ntime,scfcv_args%dtset%tolmxde) 756 end if 757 end if 758 end if 759 760 if (itime==ntime.and.icycle==ncycle) iexit=1 761 762 ! ########################################################### 763 ! ### 16. => Precondition forces, stress and energy 764 765 write(message,*) 'Geometry Optimization Precondition:',ab_mover%goprecon 766 if(need_verbose)call wrtout(std_out,message,'COLL') 767 if (ab_mover%goprecon>0)then 768 call prec_simple(ab_mover,preconforstr,hist,icycle,itime,0) 769 end if 770 771 ! ########################################################### 772 ! ### 17. => Call to each predictor 773 774 ! MT->GAF: dirty trick to predict vel(t) 775 ! do a double loop: 1- compute vel, 2- exit 776 nloop=1 777 778 779 if (scfcv_args%dtset%nctime>0.and.iexit==1) then 780 iexit=0;nloop=2 781 end if 782 783 do ii=1,nloop 784 if (ii==2) iexit=1 785 786 select case (ab_mover%ionmov) 787 case (1) 788 call pred_moldyn(ab_mover,hist,icycle,itime,ncycle,ntime,DEBUG,iexit) 789 case (2,3) 790 call pred_bfgs(ab_mover,ab_xfh,preconforstr,hist,ab_mover%ionmov,itime,DEBUG,iexit) 791 case (4,5) 792 call pred_simple(ab_mover,hist,iexit) 793 case (6,7) 794 call pred_verlet(ab_mover,hist,ab_mover%ionmov,itime,ntime,DEBUG,iexit) 795 case (8) 796 call pred_nose(ab_mover,hist,itime,ntime,DEBUG,iexit) 797 case (9) 798 call pred_langevin(ab_mover,hist,icycle,itime,ncycle,ntime,DEBUG,iexit,skipcycle) 799 case (10,11) 800 call pred_delocint(ab_mover,ab_xfh,preconforstr,hist,ab_mover%ionmov,itime,DEBUG,iexit) 801 case (12) 802 call pred_isokinetic(ab_mover,hist,itime,ntime,DEBUG,iexit) 803 case (13) 804 call pred_isothermal(ab_mover,hist,itime,mttk_vars,ntime,DEBUG,iexit) 805 case (14) 806 call pred_srkna14(ab_mover,hist,icycle,DEBUG,iexit,skipcycle) 807 case (20) 808 call pred_diisrelax(ab_mover,hist,itime,ntime,DEBUG,iexit) 809 case (21) 810 call pred_steepdesc(ab_mover,preconforstr,hist,itime,DEBUG,iexit) 811 case (22) 812 call pred_lbfgs(ab_mover,ab_xfh,preconforstr,hist,ab_mover%ionmov,itime,DEBUG,iexit) 813 #if defined HAVE_LOTF 814 case (23) 815 call pred_lotf(ab_mover,hist,itime,icycle,DEBUG,iexit) 816 #endif 817 case (24) 818 call pred_velverlet(ab_mover,hist,itime,ntime,DEBUG,iexit) 819 case (25) 820 call pred_hmc(ab_mover,hist,itime,icycle,ntime,ncycle,DEBUG,iexit) 821 822 case default 823 write(message,"(a,i0)") "Wrong value of ionmov: ",ab_mover%ionmov 824 MSG_ERROR(message) 825 end select 826 827 end do 828 829 ! check dilatmx here and correct if necessary 830 if (scfcv_args%dtset%usewvl == 0) then 831 call chkdilatmx(scfcv_args%dtset%chkdilatmx,scfcv_args%dtset%dilatmx,& 832 & rprimd,scfcv_args%dtset%rprimd_orig(1:3,1:3,1),dilatmx_errmsg) 833 _IBM6("dilatxm_errmsg: "//TRIM(dilatmx_errmsg)) 834 if (LEN_TRIM(dilatmx_errmsg) /= 0) then 835 MSG_WARNING(dilatmx_errmsg) 836 nerr_dilatmx = nerr_dilatmx+1 837 if (nerr_dilatmx > 3) then 838 ! Write last structure before aborting, so that we can restart from it. 839 ! zion is not available, but it's not useful here. 840 if (me == master) then 841 ! Init crystal 842 call crystal_init(scfcv_args%dtset%amu_orig(:,1),crystal,0,ab_mover%natom,& 843 & scfcv_args%dtset%npsp,ab_mover%ntypat,scfcv_args%dtset%nsym,rprimd,ab_mover%typat,xred,& 844 & [(-one, ii=1,ab_mover%ntypat)],ab_mover%znucl,2,.False.,.False.,"dilatmx_structure",& 845 & symrel=scfcv_args%dtset%symrel,tnons=scfcv_args%dtset%tnons,symafm=scfcv_args%dtset%symafm) 846 847 #ifdef HAVE_NETCDF 848 ! Write netcdf file 849 filename = strcat(dtfil%filnam_ds(4), "_DILATMX_STRUCT.nc") 850 NCF_CHECK(crystal_ncwrite_path(crystal, filename)) 851 #endif 852 call crystal_free(crystal) 853 end if 854 call xmpi_barrier(comm) 855 write (dilatmx_errmsg, '(a,i0,3a)') & 856 & 'Dilatmx has been exceeded too many times (', nerr_dilatmx, ')',ch10, & 857 & 'Restart your calculation from larger lattice vectors and/or a larger dilatmx' 858 MSG_ERROR_CLASS(dilatmx_errmsg, "DilatmxError") 859 end if 860 end if 861 end if 862 863 ! Write MOLDYN netcdf and POSABIN files (done every dtset%nctime time step) 864 ! This file is not created for multibinit run 865 if(need_scfcv_cycle .and. (ab_mover%ionmov/=23 .or. icycle==1))then 866 if (scfcv_args%dtset%nctime>0) then 867 jj=itime; if(hist_prev%mxhist>0.and.ab_mover%restartxf==-1) jj=jj-hist_prev%mxhist 868 if (jj>0) then 869 option=3 870 ihist_prev = abihist_findIndex(hist,-1) 871 call wrt_moldyn_netcdf(amass,scfcv_args%dtset,jj,option,dtfil%fnameabo_moldyn,& 872 & scfcv_args%mpi_enreg,scfcv_args%results_gs,& 873 & hist%rprimd(:,:,ihist_prev),dtfil%unpos,hist%vel(:,:,hist%ihist),& 874 & hist%xred(:,:,ihist_prev)) 875 end if 876 if (iexit==1) hist%ihist=ihist_prev 877 end if 878 end if 879 if(iexit/=0) exit 880 881 ! ########################################################### 882 ! ### 18. Use the history to extract the new values 883 ! ### acell, rprimd and xred 884 885 call hist2var(acell,hist,ab_mover%natom,rprimd,xred,DEBUG) 886 887 if(ab_mover%optcell/=0)then 888 889 call matr3inv(rprimd,gprimd) 890 891 ! If metric has changed since the initialization, update the Ylm's 892 if (scfcv_args%psps%useylm==1)then 893 option=0; 894 if (scfcv_args%dtset%iscf>0) option=1 895 call initylmg(gprimd,& 896 & scfcv_args%kg,& 897 & scfcv_args%dtset%kptns,& 898 & scfcv_args%dtset%mkmem,& 899 & scfcv_args%mpi_enreg,& 900 & scfcv_args%psps%mpsang,& 901 & scfcv_args%dtset%mpw,& 902 & scfcv_args%dtset%nband,& 903 & scfcv_args%dtset%nkpt,& 904 & scfcv_args%npwarr,& 905 & scfcv_args%dtset%nsppol,& 906 & option,rprimd,& 907 & scfcv_args%ylm,& 908 & scfcv_args%ylmgr) 909 end if 910 911 end if 912 913 vel(:,:)=hist%vel(:,:,hist%ihist) 914 915 ! vel_cell(3,3)= velocities of cell parameters 916 ! Not yet used here but compute it for consistency 917 vel_cell(:,:)=zero 918 if (ab_mover%ionmov==13) then 919 if (itime_hist>2) then 920 ihist_prev2 = abihist_findIndex(hist,-2) 921 vel_cell(:,:)=(hist%rprimd(:,:,hist%ihist)- hist%rprimd(:,:,ihist_prev2))/(two*ab_mover%dtion) 922 else if (itime_hist>1) then 923 ihist_prev = abihist_findIndex(hist,-1) 924 vel_cell(:,:)=(hist%rprimd(:,:,hist%ihist)-hist%rprimd(:,:,ihist_prev))/(ab_mover%dtion) 925 end if 926 end if 927 928 ! This is needed for some compilers such as 929 ! pathscale, g95, xlf that do not exit 930 ! from a loop if you change the upper limit 931 ! inside 932 if (icycle>=ncycle .and. scfcv_args%mpi_enreg%me == 0) then 933 if(need_verbose)write(std_out,*) 'EXIT:',icycle,ncycle 934 exit 935 end if 936 937 write(message,*) 'ICYCLE',icycle,skipcycle 938 if(need_verbose)call wrtout(std_out,message,'COLL') 939 write(message,*) 'NCYCLE',ncycle 940 if(need_verbose)call wrtout(std_out,message,'COLL') 941 if (skipcycle) exit 942 943 ! ########################################################### 944 ! ### 19. End loop icycle 945 946 end do ! do icycle=1,ncycle 947 948 if(iexit/=0)exit 949 950 ! ########################################################### 951 ! ### 20. End loop itime 952 953 end do ! do itime=1,ntime 954 955 ! Call fconv here if we exited due to wall time limit. 956 if (timelimit_exit==1 .and. specs%isFconv) then 957 iexit = timelimit_exit 958 ntime = itime-1 959 ihist_prev = abihist_findIndex(hist,-1) 960 if ((ab_mover%ionmov/=4.and.ab_mover%ionmov/=5)) then 961 if (scfcv_args%dtset%tolmxf/=0)then 962 call fconv(hist%fcart(:,:,ihist_prev),& 963 & scfcv_args%dtset%iatfix, & 964 & iexit, itime,& 965 & ab_mover%natom,& 966 & ntime,& 967 & ab_mover%optcell,& 968 & scfcv_args%dtset%strfact,& 969 & scfcv_args%dtset%strtarget,& 970 & hist%strten(:,ihist_prev),& 971 & scfcv_args%dtset%tolmxf) 972 else 973 call erlxconv(hist,iexit,itime,itime_hist,ntime,scfcv_args%dtset%tolmxde) 974 end if 975 end if 976 end if 977 978 ! Avoid pending requests if itime == ntime. 979 call xmpi_wait(quitsum_request,ierr) 980 981 !########################################################### 982 !### 21. Set the final values of xred with the last 983 !### computed values (not the last predicted) 984 985 call hist2var(acell,hist,ab_mover%natom,rprimd,xred,DEBUG) 986 vel(:,:)=hist%vel(:,:,hist%ihist) 987 988 if (DEBUG .and. ab_mover%ionmov==1)then 989 write (std_out,*) 'vel' 990 do kk=1,ab_mover%natom 991 write (std_out,*) hist%vel(:,kk,hist%ihist) 992 end do 993 end if 994 995 !########################################################### 996 !### 22. XML Output at the end 997 998 !XML output of the status 999 if (scfcv_args%mpi_enreg%me == 0 .and. scfcv_args%dtset%prtxml == 1) then 1000 write(ab_xml_out, "(3a)") ' <geometryMinimisation type="',trim(specs%type4xml),'">' 1001 write(ab_xml_out, "(5a)") ' <status cvState="',trim(stat4xml) ,'" stop-criterion="',trim(specs%crit4xml),'" />' 1002 write(ab_xml_out, "(3a)") ' </geometryMinimisation>' 1003 end if 1004 1005 !########################################################### 1006 !### 23. Deallocate hist and ab_mover datatypes 1007 1008 if (ab_mover%goprecon>0)then 1009 call prec_simple(ab_mover,preconforstr,hist,1,1,1) 1010 end if 1011 1012 if (ab_mover%ionmov==13)then 1013 call mttk_fin(mttk_vars) 1014 end if 1015 1016 ABI_DEALLOCATE(amu) 1017 ABI_DEALLOCATE(xred_prev) 1018 1019 call abihist_free(hist) 1020 call abihist_free(hist_prev) 1021 1022 call abimover_fin(ab_mover) 1023 call abiforstr_fin(preconforstr) 1024 1025 call status(0,dtfil%filstat,iexit,level,'exit ') 1026 1027 contains