TABLE OF CONTENTS
- ABINIT/erlxconv
- ABINIT/fconv
- ABINIT/gettag
- ABINIT/m_mover
- ABINIT/mover
- ABINIT/prtnatom
- ABINIT/prtxfase
- ABINIT/wrt_moldyn_netcdf
ABINIT/erlxconv [ Functions ]
NAME
erlxconv
FUNCTION
FIXME: add description.
INPUTS
OUTPUT
SOURCE
1166 subroutine erlxconv(hist,iexit,itime,itime_hist,ntime,tolmxde) 1167 1168 !Arguments ------------------------------------ 1169 !scalars 1170 integer,intent(in) :: itime,itime_hist,ntime 1171 integer,intent(inout) :: iexit 1172 real(dp), intent(in) :: tolmxde 1173 !arrays 1174 type(abihist),intent(inout) :: hist 1175 1176 !Local variables------------------------------- 1177 integer :: ihist,ihist_prev,ihist_prev2 1178 real(dp) :: ediff1,ediff2,maxediff 1179 character(len=500) :: msg 1180 ! ************************************************************************* 1181 1182 if (itime_hist<3) then 1183 write(msg, '(a,a,a)' ) ch10,& 1184 ' erlxconv : minimum 3 Broyd/MD steps to check convergence of energy in relaxations',ch10 1185 call wrtout(std_out,msg,'COLL') 1186 else 1187 ihist = hist%ihist 1188 ihist_prev = abihist_findIndex(hist,-1) 1189 ihist_prev2 = abihist_findIndex(hist,-2) 1190 ediff1 = hist%etot(ihist) - hist%etot(ihist_prev) 1191 ediff2 = hist%etot(ihist) - hist%etot(ihist_prev2) 1192 if ((abs(ediff1)<tolmxde).and.(abs(ediff2)<tolmxde)) then 1193 write(msg, '(a,a,i4,a,a,a,a,a,es11.4,a,a)' ) ch10,& 1194 ' At Broyd/MD step',itime,', energy is converged : ',ch10,& 1195 ' the difference in energy with respect to the two ',ch10,& 1196 ' previous steps is < tolmxde=',tolmxde,' ha',ch10 1197 call wrtout([std_out, ab_out], msg) 1198 iexit=1 1199 else 1200 maxediff = max(abs(ediff1),abs(ediff2)) 1201 if(iexit==1)then 1202 write(msg, '(a,a,a,a,i5,a,a,a,es11.4,a,es11.4,a,a)' ) ch10,& 1203 ' erlxconv : WARNING -',ch10,& 1204 ' ntime=',ntime,' was not enough Broyd/MD steps to converge energy: ',ch10,& 1205 ' max difference in energy =',maxediff,' > tolmxde=',tolmxde,' ha',ch10 1206 call wrtout([std_out, ab_out], msg) 1207 1208 write(std_out,"(8a)")ch10,& 1209 "--- !RelaxConvergenceWarning",ch10,& 1210 "message: | ",ch10,TRIM(indent(msg)),ch10,& 1211 "..." 1212 else 1213 write(msg, '(a,a,i4,a,a,a,es11.4,a,es11.4,a,a)' ) ch10,& 1214 ' erlxconv : at Broyd/MD step',itime,', energy has not converged yet. ',ch10,& 1215 ' max difference in energy=',maxediff,' > tolmxde=',tolmxde,' ha',ch10 1216 call wrtout(std_out,msg,'COLL') 1217 end if 1218 end if 1219 end if 1220 1221 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
SOURCE
1028 subroutine fconv(fcart,iatfix,iexit,itime,natom,ntime,optcell,strfact,strtarget,strten,rprim,tolmxf) 1029 1030 !Arguments ------------------------------------ 1031 !scalars 1032 integer,intent(in) :: itime,natom,ntime,optcell 1033 integer,intent(inout) :: iexit 1034 real(dp),intent(in) :: strfact,tolmxf 1035 !arrays 1036 integer,intent(in) :: iatfix(3,natom) 1037 real(dp),intent(in) :: fcart(3,natom),strtarget(6),strten(6) 1038 real(dp), intent(in) :: rprim(3,3) 1039 1040 !Local variables------------------------------- 1041 !scalars 1042 integer :: iatom,idir,istr 1043 real(dp) :: fmax,strdiag!,fcell 1044 character(len=500) :: msg 1045 !arrays 1046 real(dp) :: dstr(6) 1047 1048 ! ************************************************************************* 1049 1050 ABI_UNUSED(rprim) 1051 1052 !Compute maximal component of forces, EXCLUDING any fixed components 1053 fmax=zero 1054 do iatom=1,natom 1055 do idir=1,3 1056 if (iatfix(idir,iatom) /= 1) then 1057 if( abs(fcart(idir,iatom)) >= fmax ) fmax=abs(fcart(idir,iatom)) 1058 end if 1059 end do 1060 end do 1061 1062 dstr(:)=strten(:)-strtarget(:) 1063 1064 !Eventually take into account the stress 1065 if(optcell==1)then 1066 strdiag=(dstr(1)+dstr(2)+dstr(3))/3.0_dp 1067 if(abs(strdiag)*strfact >= fmax ) fmax=abs(strdiag)*strfact 1068 else if(optcell==2)then 1069 do istr=1,6 1070 if(abs(dstr(istr))*strfact >= fmax ) fmax=abs(dstr(istr))*strfact 1071 end do 1072 else if(optcell==3)then 1073 ! Must take away the trace from diagonal elements 1074 strdiag=(dstr(1)+dstr(2)+dstr(3))/3.0_dp 1075 do istr=1,3 1076 if(abs(dstr(istr)-strdiag)*strfact >= fmax ) fmax=abs(dstr(istr)-strdiag)*strfact 1077 end do 1078 do istr=4,6 1079 if(abs(dstr(istr))*strfact >= fmax ) fmax=abs(dstr(istr))*strfact 1080 end do 1081 ! else if(optcell==4 .or. optcell==5 .or. optcell==6)then 1082 ! if(abs(dstr(optcell-3))*strfact >= fmax ) fmax=abs(dstr(optcell-3))*strfact 1083 else if(optcell==4) then 1084 ! only dstr 1 5 6 are in xfpack_f2vout. The other component shouldn't be checked. 1085 !fcell = dstr(1) * rprim(1,1) + dstr(6) * rprim(2,1) + dstr(5) * rprim(3,1) 1086 !if (abs(fcell)*strfact >= fmax) fmax=abs(fcell)*strfact 1087 !fcell = dstr(6) * rprim(1,1) 1088 !if (abs(fcell)*strfact >= fmax) fmax=abs(fcell)*strfact 1089 !fcell = dstr(5) * rprim(1,1) 1090 !if (abs(fcell)*strfact >= fmax) fmax=abs(fcell)*strfact 1091 fmax = maxval(abs(dstr([1, 5, 6])))*strfact 1092 else if(optcell==5) then ! 2 4 6 1093 !fcell = dstr(5) * rprim(3,3) 1094 !if (abs(fcell)*strfact >= fmax) fmax=abs(fcell)*strfact 1095 !fcell = dstr(6) * rprim(1,2) + dstr(2) * rprim(2,2) + dstr(4) * rprim(3,2) 1096 !if (abs(fcell)*strfact >= fmax) fmax=abs(fcell)*strfact 1097 !fcell = dstr(4) * rprim(2,2) 1098 !if (abs(fcell)*strfact >= fmax) fmax=abs(fcell)*strfact 1099 fmax = maxval(abs(dstr([2, 4, 6])))*strfact 1100 else if(optcell==6) then 1101 !fcell = dstr(5) * rprim(3,3) 1102 !if (abs(fcell)*strfact >= fmax) fmax=abs(fcell)*strfact 1103 !fcell = dstr(4) * rprim(3,3) 1104 !if (abs(fcell)*strfact >= fmax) fmax=abs(fcell)*strfact 1105 !fcell = dstr(5) * rprim(1,3) + dstr(4) * rprim(2,3) + dstr(3) * rprim(3,3) 1106 !if (abs(fcell)*strfact >= fmax) fmax=abs(fcell)*strfact 1107 fmax = maxval(abs(dstr([3, 4, 5])))*strfact 1108 else if(optcell==7)then 1109 if(abs(dstr(2))*strfact >= fmax ) fmax=abs(dstr(2))*strfact 1110 if(abs(dstr(3))*strfact >= fmax ) fmax=abs(dstr(3))*strfact 1111 if(abs(dstr(4))*strfact >= fmax ) fmax=abs(dstr(4))*strfact 1112 else if(optcell==8)then 1113 if(abs(dstr(1))*strfact >= fmax ) fmax=abs(dstr(1))*strfact 1114 if(abs(dstr(3))*strfact >= fmax ) fmax=abs(dstr(3))*strfact 1115 if(abs(dstr(5))*strfact >= fmax ) fmax=abs(dstr(5))*strfact 1116 else if(optcell==9)then 1117 if(abs(dstr(1))*strfact >= fmax ) fmax=abs(dstr(1))*strfact 1118 if(abs(dstr(2))*strfact >= fmax ) fmax=abs(dstr(2))*strfact 1119 if(abs(dstr(6))*strfact >= fmax ) fmax=abs(dstr(6))*strfact 1120 end if 1121 1122 if (fmax<tolmxf) then 1123 write(msg, '(a,a,i4,a,a,a,es11.4,a,es11.4,a,a)' ) ch10,& 1124 ' At Broyd/MD step',itime,', gradients are converged : ',ch10,& 1125 ' max grad (force/stress) =',fmax,' < tolmxf=',tolmxf,' ha/bohr (free atoms)',ch10 1126 call wrtout([std_out, ab_out], msg) 1127 iexit=1 1128 else 1129 if(iexit==1)then 1130 write(msg, '(a,a,a,a,i5,a,a,a,es11.4,a,es11.4,a,a)' ) ch10,& 1131 ' fconv : WARNING -',ch10,& 1132 ' ntime=',ntime,' was not enough Broyd/MD steps to converge gradients: ',ch10,& 1133 ' max grad (force/stress) =',fmax,' > tolmxf=',tolmxf,' ha/bohr (free atoms)',ch10 1134 call wrtout([std_out, ab_out], msg) 1135 1136 write(std_out,"(8a)")ch10,& 1137 "--- !RelaxConvergenceWarning",ch10,& 1138 "message: | ",ch10,TRIM(indent(msg)),ch10,& 1139 "..." 1140 1141 else 1142 write(msg, '(a,i4,a,a,a,es11.4,a,es11.4,a,a)' ) & 1143 ' fconv : at Broyd/MD step',itime,', gradients have not converged yet. ',ch10,& 1144 ' max grad (force/stress) =',fmax,' > tolmxf=',tolmxf,' ha/bohr (free atoms)',ch10 1145 call wrtout(std_out,msg,'COLL') 1146 end if 1147 iexit=0 1148 end if 1149 1150 end subroutine fconv
ABINIT/gettag [ Functions ]
NAME
gettag
FUNCTION
Set the tag associated to each atom,
INPUTS
prtallatoms = Logical for PRTint ALL ATOMS atlist = ATom LIST index = index for each atom natom = Number of ATOMs
OUTPUT
tag = The string to put for each atom
SOURCE
1512 subroutine gettag(atlist,index,natom,prtallatoms,tag) 1513 1514 !Arguments ------------------------------------ 1515 !scalars 1516 logical,intent(in) :: prtallatoms 1517 integer,intent(in) :: natom 1518 logical,intent(in) :: atlist(natom) 1519 integer,intent(in) :: index 1520 character(len=7),intent(out) :: tag 1521 1522 ! ********************************************************************* 1523 !The numbering will be from (1) to (9999) 1524 1525 if (prtallatoms)then 1526 tag='' 1527 elseif (atlist(index)) then 1528 if (natom<10) then 1529 write(tag, '(a,I1.1,a)') ' (',index,')' 1530 elseif (natom<100) then 1531 write(tag, '(a,I2.2,a)') ' (',index,')' 1532 elseif (natom<1000) then 1533 write(tag, '(a,I3.3,a)') ' (',index,')' 1534 elseif (natom<10000) then 1535 write(tag, '(a,I4.4,a)') ' (',index,')' 1536 end if 1537 end if 1538 1539 end subroutine gettag
ABINIT/m_mover [ Modules ]
NAME
m_mover
FUNCTION
Move ion or change acell according to forces and stresses
COPYRIGHT
Copyright (C) 1998-2024 ABINIT group (DCA, XG, GMR, SE, FLambert,MT) 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 23 24 use defs_basis 25 use m_abicore 26 use m_errors 27 use m_profiling_abi 28 use m_abimover 29 use m_abihist 30 use m_dtset 31 use m_xmpi 32 use m_nctk 33 use m_dtfil 34 use m_yaml 35 #ifdef HAVE_NETCDF 36 use netcdf 37 #endif 38 #if defined HAVE_LOTF 39 use lotfpath 40 use m_pred_lotf 41 #endif 42 43 use defs_abitypes, only : MPI_type 44 use m_fstrings, only : strcat, sjoin, indent, itoa 45 use m_symtk, only : matr3inv, symmetrize_xred 46 use m_geometry, only : fcart2gred, chkdilatmx, xred2xcart 47 use m_time, only : abi_wtime, sec2str 48 use m_exit, only : get_start_time, have_timelimit_in, get_timelimit, enable_timelimit_in 49 use m_electronpositron, only : electronpositron_type 50 use m_scfcv, only : scfcv_t, scfcv_run 51 use m_effective_potential,only : effective_potential_type, effective_potential_evaluate 52 use m_initylmg, only : initylmg 53 use m_xfpack, only : xfh_update 54 use m_precpred_1geo, only : precpred_1geo 55 use m_pred_delocint, only : pred_delocint 56 use m_pred_bfgs, only : pred_bfgs, pred_lbfgs 57 use m_pred_fire, only : pred_fire 58 use m_pred_isokinetic, only : pred_isokinetic 59 use m_pred_diisrelax, only : pred_diisrelax 60 use m_pred_nose, only : pred_nose 61 use m_pred_srkhna14, only : pred_srkna14 62 use m_pred_isothermal, only : pred_isothermal 63 use m_pred_verlet, only : pred_verlet 64 use m_pred_velverlet, only : pred_velverlet 65 use m_pred_moldyn, only : pred_moldyn 66 use m_pred_langevin, only : pred_langevin 67 use m_pred_steepdesc, only : pred_steepdesc 68 use m_pred_simple, only : pred_simple, prec_simple 69 use m_pred_hmc, only : pred_hmc 70 !use m_generate_training_set, only : generate_training_set 71 use m_wvl_wfsinp, only : wvl_wfsinp_reformat 72 use m_wvl_rho, only : wvl_mkrho 73 use m_effective_potential_file, only : effective_potential_file_mapHistToRef 74 #if defined DEV_MS_SCALEUP 75 use scup_global, only : global_set_parent_iter,global_set_print_parameters 76 #endif 77 use m_scup_dataset 78 79 implicit none 80 81 private
ABINIT/mover [ Functions ]
NAME
mover
FUNCTION
Move ion or change acell acording to forces and stresses
INPUTS
amu_curr(ntypat)=mass of each atom for the current image 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 itimimage_gstate= [optional] counter for the itimimage loop, in the calling routine. mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol mpi_enreg=information 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) scup_dtset <type(scup_dtset_type) = derived datatype holding all options for the evaluation of an effective electronic model using SCALE UP
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 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, amu_curr, 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.
SOURCE
187 subroutine mover(scfcv_args,ab_xfh,acell,amu_curr,dtfil,& 188 & electronpositron,rhog,rhor,rprimd,vel,vel_cell,xred,xred_old,& 189 & effective_potential,filename_ddb,itimimage_gstate,verbose,writeHIST,scup_dtset,sc_size) 190 191 !Arguments ------------------------------------ 192 !scalars 193 integer, intent(in), optional :: itimimage_gstate 194 type(scfcv_t),intent(inout) :: scfcv_args 195 type(datafiles_type),intent(inout),target :: dtfil 196 type(electronpositron_type),pointer :: electronpositron 197 type(ab_xfh_type),intent(inout) :: ab_xfh 198 type(effective_potential_type),optional,intent(inout) :: effective_potential 199 logical,optional,intent(in) :: verbose 200 logical,optional,intent(in) :: writeHIST 201 character(len=fnlen),optional,intent(in) :: filename_ddb 202 !arrays 203 real(dp),intent(inout) :: acell(3) 204 real(dp), intent(in),target :: amu_curr(:) !(scfcv%dtset%ntypat) 205 real(dp), pointer :: rhog(:,:),rhor(:,:) 206 real(dp), intent(inout) :: xred(3,scfcv_args%dtset%natom),xred_old(3,scfcv_args%dtset%natom) 207 real(dp), intent(inout) :: vel(3,scfcv_args%dtset%natom),vel_cell(3,3),rprimd(3,3) 208 type(scup_dtset_type),optional, intent(inout) :: scup_dtset 209 integer,optional,intent(in) :: sc_size(3) 210 211 !Local variables------------------------------- 212 !scalars 213 integer,parameter :: level=102,master=0 214 type(abihist) :: hist,hist_prev 215 type(abimover) :: ab_mover 216 type(abimover_specs) :: specs 217 type(abiforstr) :: preconforstr ! Preconditioned forces and stress 218 type(delocint) :: deloc 219 type(mttk_type) :: mttk_vars 220 integer :: itime,icycle,itime_hist,iexit=0,ifirst,ihist_prev,ihist_prev2,timelimit_exit,ncycle,nhisttot,kk,jj,me 221 integer :: ntime,option,comm 222 integer :: nerr_dilatmx,my_quit,ierr,quitsum_request 223 integer ABI_ASYNC :: quitsum_async 224 character(len=500) :: msg 225 !character(len=500) :: dilatmx_errmsg 226 character(len=8) :: stat4xml 227 character(len=35) :: fmt 228 character(len=fnlen) :: filename,fname_ddb,name_file 229 character(len=500) :: MY_NAME = "mover" 230 real(dp) :: gr_avg 231 logical :: DEBUG=.FALSE., need_verbose=.TRUE.,need_writeHIST=.TRUE. 232 logical :: need_scfcv_cycle = .TRUE., need_elec_eval = .FALSE. 233 logical :: changed,useprtxfase 234 logical :: skipcycle 235 integer :: minIndex,ii,similar,conv_retcode 236 integer :: iapp 237 logical :: file_exists 238 real(dp) :: minE,wtime_step,now,prev 239 !arrays 240 integer :: itimes(2) 241 real(dp) :: gprimd(3,3),rprim(3,3),rprimd_prev(3,3) 242 real(dp),allocatable :: gred_corrected(:,:),xred_prev(:,:) 243 ! *************************************************************** 244 need_verbose=.TRUE. 245 if(present(verbose)) need_verbose = verbose 246 247 need_writeHIST=.TRUE. 248 if(present(writeHIST)) need_writeHIST = writeHIST 249 250 ! enable time limit handler if not done in callers. 251 if (enable_timelimit_in(MY_NAME) == MY_NAME) then 252 if (need_verbose) then 253 write(std_out,*)"Enabling timelimit check in function: ",trim(MY_NAME)," with timelimit: ",trim(sec2str(get_timelimit())) 254 end if 255 end if 256 257 !Table of contents 258 !(=>) Refers to an important call (scfcv,pred_*) 259 ! 260 !01. Initialization of indexes and allocations of arrays 261 !02. Particularities of each predictor 262 !03. Set the number of iterations ntime 263 !04. Try to read history of previous calculations 264 !05. Allocate the hist structure 265 !06. First output before any itime or icycle 266 !07. Fill the history of the first SCFCV 267 !08. Loop for itime (From 1 to ntime) 268 !09. Loop for icycle (From 1 to ncycle) 269 !10. Output for each icycle (and itime) 270 !11. Symmetrize atomic coordinates over space group elements 271 !12. => Call to SCFCV routine and fill history with forces 272 !13. Write the history into the _HIST file 273 !14. Output after SCFCV 274 !15. => Test Convergence of forces and stresses 275 !16. => Precondition forces, stress and energy 276 !17. => Call to each predictor 277 !18. Use the history to extract the new values 278 !19. End loop icycle 279 !20. End loop itime 280 !21. Set the final values of xred 281 !22. XML Output at the end 282 !23. Deallocate hist and ab_mover datatypes 283 ! 284 call abimover_ini(ab_mover,amu_curr,dtfil,scfcv_args%dtset,specs) 285 286 if (ab_mover%ionmov==10 .or. ab_mover%ionmov==11) call delocint_ini(deloc) 287 288 if (ab_mover%ionmov==13 .or. ab_mover%ionmov==25)then 289 call mttk_ini(mttk_vars,ab_mover%nnos) 290 end if 291 292 !########################################################### 293 !### 03. Set the number of iterations ntime 294 !### By default ntime==1 but if the user enters a lower 295 !### value, mover will execute at least one iteration 296 297 if (scfcv_args%dtset%ntime<1)then 298 ntime=1 299 else 300 ntime=scfcv_args%dtset%ntime 301 end if 302 303 !########################################################### 304 !### 04. Try to read history of previous calculations 305 !### It requires NetCDF library 306 307 !Init MPI data 308 comm=scfcv_args%mpi_enreg%comm_cell 309 me=xmpi_comm_rank(comm) 310 311 312 #if defined HAVE_NETCDF 313 filename=trim(ab_mover%filnam_ds(4))//'_HIST.nc' 314 315 if (ab_mover%restartxf<0)then 316 ! Read history from file (and broadcast if MPI) 317 if (me==master) then 318 call read_md_hist(filename,hist_prev,specs%isVused,specs%isARused,ab_mover%restartxf==-3) 319 end if 320 call abihist_bcast(hist_prev,master,comm) 321 322 ! If restartxf specifies to reconstruct the history 323 if (hist_prev%mxhist>0.and.ab_mover%restartxf==-1)then 324 ntime=ntime+hist_prev%mxhist 325 end if 326 327 ! If restartxf specifies to start from the lowest energy 328 if (hist_prev%mxhist>0.and.ab_mover%restartxf==-2)then 329 minE=hist_prev%etot(1) 330 minIndex=1 331 do ii=1,hist_prev%mxhist 332 if(need_verbose) write(std_out,*) 'Iteration:',ii,' Total Energy:',hist_prev%etot(ii) 333 if (minE>hist_prev%etot(ii))then 334 minE=hist_prev%etot(ii) 335 minIndex=ii 336 end if 337 end do 338 if(need_verbose)write(std_out,*) 'The lowest energy occurs at iteration:',minIndex,'etotal=',minE 339 acell(:) =hist_prev%acell(:,minIndex) 340 rprimd(:,:)=hist_prev%rprimd(:,:,minIndex) 341 xred(:,:) =hist_prev%xred(:,:,minIndex) 342 vel(:, :) = hist_prev%vel(:, :, minIndex) 343 call abihist_free(hist_prev) 344 end if 345 ! If restarxf specifies to start to the last iteration 346 if (hist_prev%mxhist>0.and.ab_mover%restartxf==-3)then 347 if(present(effective_potential))then 348 call effective_potential_file_mapHistToRef(effective_potential,hist_prev,comm,scfcv_args%dtset%iatfix,need_verbose,sc_size) ! Map Hist to Ref to order atoms 349 !xred(:,:) = hist_prev%xred(:,:,1) ! Fill xred with new ordering 350 hist%ihist = 1 351 end if 352 acell(:) =hist_prev%acell(:,hist_prev%mxhist) 353 rprimd(:,:)=hist_prev%rprimd(:,:,hist_prev%mxhist) 354 xred(:,:) =hist_prev%xred(:,:,hist_prev%mxhist) 355 vel(:, :) = hist_prev%vel(:, :, hist_prev%mxhist) 356 call abihist_free(hist_prev) 357 end if 358 359 end if !if (ab_mover%restartxf<=0) 360 #endif 361 362 !########################################################### 363 !### 05. Allocate the hist structure 364 365 iexit=0; timelimit_exit=0 366 ncycle=specs%ncycle 367 368 if(ab_mover%ionmov==25.and.scfcv_args%dtset%hmctt>=0)then 369 ncycle=scfcv_args%dtset%hmctt 370 if(scfcv_args%dtset%hmcsst>0.and.ab_mover%optcell/=0)then 371 ncycle=ncycle+scfcv_args%dtset%hmcsst 372 endif 373 endif 374 375 nhisttot=ncycle*ntime;if (scfcv_args%dtset%nctime>0) nhisttot=nhisttot+1 376 !AM_2017 New version of the hist, we just store the needed history step not all of them... 377 if(specs%nhist/=-1)then 378 nhisttot = specs%nhist! We don't need to store all the history 379 endif 380 381 call abihist_init(hist,ab_mover%natom,nhisttot,specs%isVused,specs%isARused) 382 call abiforstr_ini(preconforstr,ab_mover%natom) 383 384 !########################################################### 385 !### 06. First output before any itime or icycle iteration 386 387 !If effective potential is present forces will be compute with it 388 if (present(effective_potential))then 389 need_scfcv_cycle = .FALSE. 390 if(need_verbose)then 391 write(msg,'(2a,i2,5a,80a)')& 392 & ch10,'=== [ionmov=',ab_mover%ionmov,'] ',trim(specs%method),' with effective potential',& 393 & ch10,('=',kk=1,80) 394 call wrtout([std_out, ab_out], msg) 395 end if 396 need_elec_eval = .FALSE. 397 if(present(scup_dtset))then 398 need_elec_eval = scup_dtset%scup_elec_model 399 endif 400 else 401 if(need_verbose)then 402 write(msg,'(a,a,i2,a,a,a,80a)')& 403 & ch10,'=== [ionmov=',ab_mover%ionmov,'] ',specs%method,& 404 & ch10,('=',kk=1,80) 405 call wrtout([std_out, ab_out], msg) 406 end if 407 end if 408 409 !Format for printing on each cycle 410 write(fmt,'(a6,i2,a4,i2,a4,i2,a4,i2,a9)')& 411 '(a,a,i',int(log10(real(ntime))+1),& 412 ',a,i',int(log10(real(ntime))+1),& 413 ',a,i',int(log10(real(ncycle))+1),& 414 ',a,i',int(log10(real(ncycle))+1),& 415 ',a,a,80a)' 416 417 !########################################################### 418 !### 07. Fill the history of the first SCFCV 419 420 if (ab_mover%ionmov==26)then 421 422 !Tdep call need to merge with adewandre branch 423 else if (ab_mover%ionmov==27)then 424 if(present(filename_ddb))then 425 fname_ddb = trim(filename_ddb) 426 else 427 fname_ddb = trim(ab_mover%filnam_ds(3))//'_DDB' 428 end if 429 INQUIRE(FILE=filename, EXIST=file_exists) 430 431 ABI_ERROR("This section has been disabled, ph_freez_disp is not defined in main ABINIT") 432 433 ! XG 20200322 : The input variables ph_freez_disp are not documented neither tested, so they 434 ! have been removed from the allowed list in the parser. Also, you should not be here ! 435 ! call generate_training_set(acell,ab_mover%ph_freez_disp_addStrain==1,ab_mover%ph_freez_disp_ampl,& 436 !& fname_ddb,hist,ab_mover%natom,ab_mover%ph_freez_disp_nampl,ntime,& 437 !& ab_mover%ph_ngqpt,ab_mover%ph_nqshift,ab_mover%ph_freez_disp_option,& 438 !& ab_mover%ph_qshift,scfcv_args%dtset%supercell_latt,& 439 !& rprimd,ab_mover%mdtemp(2),xred,comm,DEBUG) 440 441 442 !Fill history with the values of xred, acell and rprimd of the first configuration 443 acell(:) =hist%acell(:,1) 444 rprimd(:,:)=hist%rprimd(:,:,1) 445 xred(:,:) =hist%xred(:,:,1) 446 447 else 448 449 ! Fill history with the values of xred, acell and rprimd 450 call var2hist(acell,hist,ab_mover%natom,rprimd,xred,DEBUG) 451 452 ! Fill velocities and ionic kinetic energy 453 call vel2hist(ab_mover%amass,hist,vel,vel_cell) 454 hist%time(hist%ihist)=zero 455 456 end if 457 458 !Decide if prtxfase will be called 459 useprtxfase=.FALSE. 460 do ii=1,ab_mover%natom 461 if (ab_mover%prtatlist(ii)/=0)then 462 useprtxfase=.TRUE.; exit 463 end if 464 end do 465 466 !At beginning no error 467 nerr_dilatmx = 0 468 469 ABI_MALLOC(xred_prev,(3,scfcv_args%dtset%natom)) 470 471 !########################################################### 472 !### 08. Loop for itime (From 1 to ntime) 473 quitsum_request = xmpi_request_null 474 475 do itime=1,ntime 476 477 call yaml_iterstart("itime", itime, dev_null, scfcv_args%dtset%use_yaml) 478 479 ! Handle time limit condition. 480 if (itime == 1) prev = abi_wtime() 481 if (itime > 1) then 482 now = abi_wtime() 483 wtime_step = now - prev 484 prev = now 485 write(msg,*)sjoin("{mover_itime:", itoa(itime - 1), ", wall_time: '", sec2str(wtime_step), "'} <<< TIME") 486 if(need_verbose)call wrtout(std_out, msg) 487 if (have_timelimit_in(MY_NAME)) then 488 if (itime > 2) then 489 call xmpi_wait(quitsum_request,ierr) 490 if (quitsum_async > 0) then 491 write(msg,"(3a)")"Approaching time limit ",trim(sec2str(get_timelimit())), ". Will exit itime loop in mover." 492 if(need_verbose) ABI_COMMENT(msg) 493 if(need_verbose) call wrtout(ab_out, msg) 494 timelimit_exit = 1 495 exit 496 end if 497 end if 498 499 my_quit = 0; if (now - get_start_time() + 2.15 * wtime_step > get_timelimit()) my_quit = 1 500 call xmpi_isum(my_quit,quitsum_async,comm,quitsum_request,ierr) 501 end if 502 end if 503 504 skipcycle=.FALSE. 505 #if defined HAVE_LOTF 506 if(ab_mover%ionmov==23 .and. .not. lotf_extrapolation(itime)) skipcycle=.True. 507 #endif 508 509 ! If RMM-DIIS is used, decrease the number of NSCF steps done with wfoptalg before activating RMM-DIIS. 510 ! In vtowfk we have the condition: istep > 3 + dtset%rmm_diis 511 ! so setting rmm_diis = 1 gives: 512 ! 4 NSCF iterations for itime == 1 513 ! 1 NSCF iterations for itime >= 2. 514 if (scfcv_args%dtset%rmm_diis /= 0 .and. itime == 2) then 515 scfcv_args%dtset%rmm_diis = scfcv_args%dtset%rmm_diis - 3 516 if (scfcv_args%dtset%rmm_diis == 0) scfcv_args%dtset%rmm_diis = 1 517 call wrtout(std_out, sjoin(" itime == 2 with RMM-DIIS --> setting rmm_diis to:", itoa(scfcv_args%dtset%rmm_diis))) 518 end if 519 520 ! ########################################################### 521 ! ### 09. Loop for icycle (From 1 to ncycle) 522 do icycle=1,ncycle 523 524 call yaml_iterstart("icycle", icycle, dev_null, scfcv_args%dtset%use_yaml) 525 526 itime_hist = (itime-1)*ncycle + icycle ! Store the time step in the history 527 528 ! ########################################################### 529 ! ### 10. Output for each icycle (and itime) 530 if(need_verbose)then 531 write(msg,fmt)& 532 ch10,'--- Iteration: (',itime,'/',ntime,') Internal Cycle: (',icycle,'/',ncycle,')',ch10,('-',kk=1,80) 533 call wrtout([std_out, ab_out], msg) 534 end if 535 if (useprtxfase) call prtxfase(ab_mover,hist,itime_hist,std_out,mover_BEFORE) 536 537 xred_prev(:,:)=xred(:,:) 538 rprimd_prev(:,:)=rprimd(:,:) 539 540 ! ########################################################### 541 ! ### 11. Symmetrize atomic coordinates over space group elements 542 543 call symmetrize_xred(ab_mover%natom,& 544 scfcv_args%dtset%nsym,scfcv_args%dtset%symrel,scfcv_args%dtset%tnons,xred,indsym=scfcv_args%indsym) 545 546 changed = any(xred /= xred_prev) 547 if (changed)then 548 hist%xred(:,:,hist%ihist)=xred(:,:) 549 if(need_verbose) then 550 write(std_out,*) 'WARNING: ATOMIC COORDINATES WERE SYMMETRIZED' 551 write(std_out,*) 'DIFFERENCES:' 552 do kk=1,ab_mover%natom 553 write(std_out,*) xred(:,kk)-xred_prev(:,kk) 554 end do 555 end if 556 xred_prev(:,:)=xred(:,:) 557 end if 558 559 ! ########################################################### 560 ! ### 12. => Call to SCFCV routine and fill history with forces 561 if (need_verbose) then 562 if (need_scfcv_cycle) then 563 write(msg,'(a,3a,33a,44a)')& 564 ch10,('-',kk=1,3),'SELF-CONSISTENT-FIELD CONVERGENCE',('-',kk=1,44) 565 else 566 write(msg,'(a,3a,33a,44a)')& 567 ch10,('-',kk=1,3),'EFFECTIVE POTENTIAL CALCULATION',('-',kk=1,44) 568 end if 569 call wrtout([std_out, ab_out], msg) 570 end if 571 572 if(hist_prev%mxhist>0.and.ab_mover%restartxf==-1.and.hist_prev%ihist<=hist_prev%mxhist)then 573 574 call abihist_compare_and_copy(hist_prev,hist,ab_mover%natom,similar,tol8,specs%nhist==nhisttot) 575 hist_prev%ihist=hist_prev%ihist+1 576 577 else 578 scfcv_args%ndtpawuj=0 579 iapp=itime 580 if(icycle>1.and.icycle/=ncycle) iapp=-1 581 if(itime==1 .and. icycle/=ncycle ) iapp=-icycle-1 582 if (ab_mover%ionmov==14.and.(icycle<ncycle)) iapp=-1 583 584 #if defined HAVE_LOTF 585 if (ab_mover%ionmov/=23 .or.(lotf_extrapolation(itime).and.(icycle/=1.or.itime==1)))then 586 #endif 587 !call scfcv_new2(scfcv_args,electronpositron,rhog,rhor,rprimd,xred,xred_old,conv_retcode) 588 589 !WVL - reformat the wavefunctions in the case of xred != xred_old 590 if (scfcv_args%dtset%usewvl == 1 .and. maxval(xred_old - xred) > zero) then 591 ! Before running scfcv, on non-first geometry step iterations, 592 ! we need to reformat the wavefunctions, taking into acount the new 593 ! coordinates. We prepare to change rhog (to be removed) and rhor. 594 ABI_FREE(rhog) 595 ABI_FREE(rhor) 596 call wvl_wfsinp_reformat(scfcv_args%dtset, scfcv_args%mpi_enreg,& 597 & scfcv_args%psps, rprimd, scfcv_args%wvl, xred, xred_old) 598 scfcv_args%nfftf = scfcv_args%dtset%nfft 599 ABI_MALLOC(rhog,(2, scfcv_args%dtset%nfft)) 600 ABI_MALLOC(rhor,(2, scfcv_args%dtset%nfft)) 601 call wvl_mkrho(scfcv_args%dtset, scfcv_args%irrzon, scfcv_args%mpi_enreg,& 602 & scfcv_args%phnons, rhor,scfcv_args%wvl%wfs,scfcv_args%wvl%den) 603 end if 604 605 ! MAIN CALL TO SELF-CONSISTENT FIELD ROUTINE 606 if (need_scfcv_cycle) then 607 608 call dtfil_init_time(dtfil,iapp) 609 itimes(1)=itime ; itimes(2)=1 610 if(present(itimimage_gstate))then 611 itimes(2)=itimimage_gstate 612 endif 613 614 !DEBUG 615 ! write(std_out,'(a,5i4)')' m_mover, before scfcv_run : itimes(1:2)=',itimes(1:2) 616 !ENDDEBUG 617 call scfcv_run(scfcv_args, electronpositron, itimes, rhog, rhor, rprimd, xred, xred_old, conv_retcode) 618 if (conv_retcode == -1) then 619 msg = "Scf cycle returned conv_retcode == -1 (timelimit is approaching), this should not happen inside mover" 620 ABI_WARNING(msg) 621 end if 622 623 else 624 ! For monte carlo don't need to recompute energy here (done in pred_montecarlo) 625 name_file='MD_anharmonic_terms_energy.dat' 626 if(itime == 1 .and. ab_mover%restartxf==-3)then 627 if(icycle==1)call effective_potential_file_mapHistToRef(effective_potential,hist,comm,scfcv_args%dtset%iatfix,& 628 & need_verbose,sc_size=sc_size)!Map Hist to Ref to order atoms 629 xred(:,:) = hist%xred(:,:,1) ! Fill xred with new ordering 630 hist%ihist = 1 631 end if 632 633 #if defined DEV_MS_SCALEUP 634 !If we a SCALE UP effective electron model give the iteration and set print-options 635 if(need_elec_eval)then 636 call global_set_parent_iter(itime) 637 ! Set all print options to false. 638 call global_set_print_parameters(geom=.FALSE.,eigvals=.FALSE.,eltic=.FALSE.,& 639 & orbocc=.FALSE.,bands=.FALSE.) 640 if(itime == 1 .or. modulo(itime,scup_dtset%scup_printniter) == 0)then 641 call global_set_print_parameters(scup_dtset%scup_printgeom,scup_dtset%scup_printeigv,scup_dtset%scup_printeltic,& 642 & scup_dtset%scup_printorbocc,scup_dtset%scup_printbands) 643 end if 644 end if 645 #endif 646 647 call effective_potential_evaluate( & 648 & effective_potential,scfcv_args%results_gs%etotal,scfcv_args%results_gs%fcart,scfcv_args%results_gs%gred,& 649 & scfcv_args%results_gs%strten,ab_mover%natom,rprimd,xred=xred,verbose=need_verbose,& 650 & filename=name_file,elec_eval=need_elec_eval) 651 652 ! Check if the simulation did not diverge... 653 if(itime > 3 .and.ABS(scfcv_args%results_gs%etotal - hist%etot(1)) > 1E5)then 654 ! We set to false the flag corresponding to the bound 655 effective_potential%anharmonics_terms%bounded = .FALSE. 656 if(need_verbose.and.me==master)then 657 ABI_WARNING("The simulation is diverging, please check your effective potential") 658 end if 659 ! Set the flag to finish the simulation 660 iexit=1 661 stat4xml="Failed" 662 else 663 ! We set to true the flag corresponding to the bound 664 effective_potential%anharmonics_terms%bounded = .TRUE. 665 end if 666 end if 667 #if defined HAVE_LOTF 668 end if 669 #endif 670 ! ANOMALOUS SITUATION 671 ! This is the only case where rprimd could change inside scfcv 672 ! It generates a weird condition, we start with a certain 673 ! value for rprimd before scfcv and after we finish with a different value. 674 ! Notice that normally scfcv should not change rprimd 675 ! And even worse if optcell==0 676 ! The solution here is to recompute acell and store these value 677 ! in the present record even if initially it was not exactly 678 ! the value entering in scfcv 679 ! One test case with these condition is bigdft/t10 680 if (any(rprimd(:,:)/=rprimd_prev(:,:))) then 681 hist%acell(:,hist%ihist)=acell(:) 682 hist%rprimd(:,:,hist%ihist)=rprimd(:,:) 683 end if 684 685 ! ANOMALOUS SITUATIONS 686 ! * In ionmov 4 & 5 xred could change inside SCFCV 687 ! So we need to take the values from the output 688 ! 689 ! * Inside scfcv_core.F90 there is a call to symmetrize_xred.F90 690 ! for the first SCF cycle symmetrize_xred could change xred 691 if (ab_mover%ionmov<10)then 692 changed = any(xred /= xred_prev) 693 if (changed)then 694 hist%xred(:,:,hist%ihist)=xred(:,:) 695 if(need_verbose)then 696 write(std_out,*) 'WARNING: ATOMIC COORDINATES WERE SYMMETRIZED AFTER SCFCV' 697 write(std_out,*) 'DIFFERENCES:' 698 do kk=1,ab_mover%natom 699 write(std_out,*) xred(:,kk)-xred_prev(:,kk) 700 end do 701 end if 702 end if 703 end if 704 705 ! Fill velocities and ionic kinetic energy 706 call vel2hist(ab_mover%amass,hist,vel,vel_cell) 707 708 hist%fcart(:,:,hist%ihist)=scfcv_args%results_gs%fcart(:,:) 709 hist%strten(:,hist%ihist) =scfcv_args%results_gs%strten(:) 710 hist%etot(hist%ihist) =scfcv_args%results_gs%etotal 711 hist%entropy(hist%ihist) =scfcv_args%results_gs%energies%entropy 712 hist%time(hist%ihist) =real(itime,kind=dp) 713 714 ! !###################################################################### 715 end if ! if(hist_prev%mxhist>0.and.ab_mover%restartxf==-1.and.hist_prev%ihist<=hist_prev%mxhist)then 716 717 ! Store trajectory in xfh file 718 if((ab_xfh%nxfh==0.or.itime/=1)) then 719 ABI_MALLOC(gred_corrected,(3,scfcv_args%dtset%natom)) 720 call fcart2gred(hist%fcart(:,:,hist%ihist),gred_corrected,rprimd,ab_mover%natom) 721 ! Get rid of mean force on whole unit cell, 722 ! but only if no generalized constraints are in effect 723 if (ab_mover%nconeq==0)then 724 do ii=1,3 725 if (ii/=3.or.ab_mover%jellslab==0) then 726 gr_avg=sum(gred_corrected(ii,:))/dble(ab_mover%natom) 727 gred_corrected(ii,:)=gred_corrected(ii,:)-gr_avg 728 end if 729 end do 730 end if 731 if (ncycle<10.and.ab_mover%restartxf>=0) then 732 do ii=1,3 733 rprim(ii,1:3)=rprimd(ii,1:3)/acell(1:3) 734 end do 735 736 ! The size of ab_xfh%xfhist is to big for very large supercell. 737 ! Call it only for specific ionmov 738 if(any((/2,3,10,11,22/)==ab_mover%ionmov)) then 739 call xfh_update(ab_xfh,acell,gred_corrected,ab_mover%natom,rprim,hist%strten(:,hist%ihist),xred) 740 end if 741 end if 742 ABI_FREE(gred_corrected) 743 end if 744 745 ! ########################################################### 746 ! ### 13. Write the history into the _HIST file 747 ! ### 748 749 #if defined HAVE_NETCDF 750 if (need_writeHIST.and.me==master) then 751 ifirst=merge(0,1,(itime>1.or.icycle>1)) 752 call write_md_hist(hist,filename,ifirst,itime_hist,ab_mover%natom,scfcv_args%dtset%nctime,& 753 & ab_mover%ntypat,ab_mover%typat,amu_curr,ab_mover%znucl,ab_mover%dtion,scfcv_args%dtset%mdtemp) 754 end if 755 #endif 756 757 ! ########################################################### 758 ! ### 14. Output after SCFCV 759 if(need_verbose.and.need_scfcv_cycle)then 760 write(msg,'(a,3a,a,72a)')ch10,('-',kk=1,3),'OUTPUT',('-',kk=1,71) 761 call wrtout([std_out, ab_out], msg) 762 end if 763 if (useprtxfase) then 764 call prtxfase(ab_mover,hist,itime_hist,ab_out,mover_AFTER) 765 call prtxfase(ab_mover,hist,itime_hist,std_out,mover_AFTER) 766 end if 767 768 ! ########################################################### 769 ! ### 15. => Test Convergence of forces and stresses 770 771 if (itime==ntime.and.icycle==ncycle)then 772 iexit=1 773 stat4xml="Failed" 774 else 775 stat4xml="Succeded" 776 end if 777 778 ! Only if convergence is needed 779 if(specs%isFconv)then 780 if ((ab_mover%ionmov/=4.and.ab_mover%ionmov/=5).or.mod(itime,2)==1)then 781 if (scfcv_args%dtset%tolmxf/=0)then 782 783 call fconv(hist%fcart(:,:,hist%ihist),& 784 & scfcv_args%dtset%iatfix, & 785 & iexit,itime,& 786 & ab_mover%natom,& 787 & ntime,& 788 & ab_mover%optcell,& 789 & scfcv_args%dtset%strfact,& 790 & scfcv_args%dtset%strtarget,& 791 & hist%strten(:,hist%ihist),& 792 & rprim,& 793 & scfcv_args%dtset%tolmxf) 794 else 795 call erlxconv(hist,iexit,itime,itime_hist,ntime,scfcv_args%dtset%tolmxde) 796 end if 797 end if 798 end if 799 800 if (itime==ntime.and.icycle==ncycle) iexit=1 801 802 ! ########################################################### 803 ! ### 16. => Precondition forces, stress and energy 804 ! ### 17. => Call to each predictor 805 806 call precpred_1geo(ab_mover,ab_xfh,amu_curr,deloc,& 807 & scfcv_args%dtset%chkdilatmx,& 808 & scfcv_args%mpi_enreg%comm_cell,& 809 & scfcv_args%dtset%dilatmx,dtfil%filnam_ds(4),& 810 & hist,scfcv_args%dtset%hmctt,& 811 & icycle,iexit,itime,mttk_vars,& 812 & scfcv_args%dtset%nctime,ncycle,nerr_dilatmx,scfcv_args%dtset%npsp,ntime,& 813 & scfcv_args%dtset%rprimd_orig,skipcycle,& 814 & scfcv_args%dtset%usewvl) 815 816 ! Write MOLDYN netcdf and POSABIN files (done every dtset%nctime time step) 817 ! This file is not created for multibinit run 818 if(need_scfcv_cycle .and. (ab_mover%ionmov/=23 .or. icycle==1))then 819 if (scfcv_args%dtset%nctime>0) then 820 jj=itime; if(hist_prev%mxhist>0.and.ab_mover%restartxf==-1) jj=jj-hist_prev%mxhist 821 if (jj>0) then 822 option=3 823 ihist_prev = abihist_findIndex(hist,-1) 824 call wrt_moldyn_netcdf(ab_mover%amass,scfcv_args%dtset,jj,option,dtfil%fnameabo_moldyn,& 825 & scfcv_args%mpi_enreg,scfcv_args%results_gs,& 826 & hist%rprimd(:,:,ihist_prev),dtfil%unpos,hist%vel(:,:,hist%ihist),& 827 & hist%xred(:,:,ihist_prev)) 828 end if 829 if (iexit==1) hist%ihist=ihist_prev 830 end if 831 end if 832 if(iexit/=0) exit 833 834 ! ########################################################### 835 ! ### 18. Use the history to extract the new values of acell, rprimd and xred 836 837 call hist2var(acell,hist,ab_mover%natom,rprimd,xred,DEBUG) 838 839 if (ab_mover%optcell/=0) then 840 ! Cell may change 841 842 call matr3inv(rprimd,gprimd) 843 844 ! If metric has changed since the initialization, update the Ylm's 845 if (scfcv_args%psps%useylm==1)then 846 option=0 847 if (scfcv_args%dtset%iscf>0) option=1 848 call initylmg(gprimd,& 849 & scfcv_args%kg,& 850 & scfcv_args%dtset%kptns,& 851 & scfcv_args%dtset%mkmem,& 852 & scfcv_args%mpi_enreg,& 853 & scfcv_args%psps%mpsang,& 854 & scfcv_args%dtset%mpw,& 855 & scfcv_args%dtset%nband,& 856 & scfcv_args%dtset%nkpt,& 857 & scfcv_args%npwarr,& 858 & scfcv_args%dtset%nsppol,& 859 & option,rprimd,& 860 & scfcv_args%ylm,& 861 & scfcv_args%ylmgr) 862 end if 863 864 end if 865 866 vel(:,:)=hist%vel(:,:,hist%ihist) 867 868 ! vel_cell(3,3)= velocities of cell parameters 869 ! Not yet used here but compute it for consistency 870 vel_cell(:,:)=zero 871 if (ab_mover%ionmov==13 .and. hist%mxhist >= 2) then 872 if (itime_hist>2) then 873 ihist_prev2 = abihist_findIndex(hist,-2) 874 vel_cell(:,:)=(hist%rprimd(:,:,hist%ihist)- hist%rprimd(:,:,ihist_prev2))/(two*ab_mover%dtion) 875 else if (itime_hist>1) then 876 ihist_prev = abihist_findIndex(hist,-1) 877 vel_cell(:,:)=(hist%rprimd(:,:,hist%ihist)-hist%rprimd(:,:,ihist_prev))/(ab_mover%dtion) 878 end if 879 end if 880 881 ! This is needed for some compilers such as 882 ! pathscale, g95, xlf that do not exit 883 ! from a loop if you change the upper limit 884 ! inside 885 if (icycle>=ncycle .and. scfcv_args%mpi_enreg%me == 0) then 886 if(need_verbose)write(std_out,*) 'EXIT:',icycle,ncycle 887 exit 888 end if 889 890 891 if (need_verbose) then 892 write(msg,*) 'ICYCLE',icycle,skipcycle 893 call wrtout(std_out,msg) 894 write(msg,*) 'NCYCLE',ncycle 895 call wrtout(std_out,msg) 896 end if 897 if (skipcycle) exit 898 899 ! ########################################################### 900 ! ### 19. End loop icycle 901 902 end do ! icycle 903 904 if(iexit/=0)exit 905 906 ! ########################################################### 907 ! ### 20. End loop itime 908 909 end do ! itime 910 911 ! Call fconv here if we exited due to wall time limit. 912 if (timelimit_exit==1 .and. specs%isFconv) then 913 iexit = timelimit_exit 914 ntime = itime-1 915 ihist_prev = abihist_findIndex(hist,-1) 916 if ((ab_mover%ionmov/=4.and.ab_mover%ionmov/=5)) then 917 if (scfcv_args%dtset%tolmxf/=0)then 918 call fconv(hist%fcart(:,:,ihist_prev),& 919 & scfcv_args%dtset%iatfix, & 920 & iexit, itime,& 921 & ab_mover%natom,& 922 & ntime,& 923 & ab_mover%optcell,& 924 & scfcv_args%dtset%strfact,& 925 & scfcv_args%dtset%strtarget,& 926 & hist%strten(:,ihist_prev),& 927 & rprim,& 928 & scfcv_args%dtset%tolmxf) 929 else 930 call erlxconv(hist,iexit,itime,itime_hist,ntime,scfcv_args%dtset%tolmxde) 931 end if 932 end if 933 end if 934 935 ! Avoid pending requests if itime == ntime. 936 call xmpi_wait(quitsum_request,ierr) 937 938 !########################################################### 939 !### 21. Set the final values of xred with the last computed values (not the last predicted) 940 941 call hist2var(acell,hist,ab_mover%natom,rprimd,xred,DEBUG) 942 vel(:,:)=hist%vel(:,:,hist%ihist) 943 944 if (DEBUG .and. ab_mover%ionmov==1)then 945 write (std_out,*) 'vel' 946 do kk=1,ab_mover%natom 947 write (std_out,*) hist%vel(:,kk,hist%ihist) 948 end do 949 end if 950 951 !########################################################### 952 !### 22. XML Output at the end 953 954 !XML output of the status 955 if (scfcv_args%mpi_enreg%me == 0 .and. scfcv_args%dtset%prtxml == 1) then 956 write(ab_xml_out, "(3a)") ' <geometryMinimisation type="',trim(specs%type4xml),'">' 957 write(ab_xml_out, "(5a)") ' <status cvState="',trim(stat4xml) ,'" stop-criterion="',trim(specs%crit4xml),'" />' 958 write(ab_xml_out, "(3a)") ' </geometryMinimisation>' 959 end if 960 961 !########################################################### 962 !### 23. Deallocate hist and ab_mover datatypes 963 964 !This call is needed to free an internal matrix. However, this is not optimal ... 965 !One should instead have a datastructure associated with the preconditioner. 966 if (ab_mover%goprecon>0) call prec_simple(ab_mover,preconforstr,hist,1,1,1) 967 968 if (ab_mover%ionmov==13 .or. ab_mover%ionmov==25)then 969 call mttk_fin(mttk_vars) 970 end if 971 972 if (ab_mover%ionmov==10 .or. ab_mover%ionmov==11) call delocint_fin(deloc) 973 974 ABI_FREE(xred_prev) 975 976 call abihist_free(hist) 977 call abihist_free(hist_prev) 978 call abimover_destroy(ab_mover) 979 call abiforstr_fin(preconforstr) 980 981 contains
ABINIT/prtnatom [ Functions ]
NAME
prtnatom
FUNCTION
Print information for N atoms
INPUTS
prtallatoms = Logical for PRTint ALL ATOMS atlist = ATom LIST index = index for each atom natom = Number of ATOMs
OUTPUT
tag = The string to put for aech atom
SOURCE
1562 subroutine prtnatom(atlist,iout,message,natom,prtallatoms,thearray) 1563 1564 !Arguments ------------------------------------ 1565 !scalars 1566 logical,intent(in) :: prtallatoms 1567 integer,intent(in) :: natom 1568 logical,intent(in) :: atlist(natom) 1569 integer,intent(in) :: iout 1570 character(len=*),intent(inout) :: message 1571 !arrays 1572 real(dp) :: thearray(3,natom) 1573 1574 !Local variables------------------------------- 1575 !scalars 1576 integer :: kk 1577 character(len=7) :: tag ! Maximal ' (9999)' 1578 character(len=18) :: fmt 1579 1580 ! ********************************************************************* 1581 1582 fmt='(a,a,1p,3e22.14,a)' 1583 1584 do kk=1,natom 1585 if (atlist(kk)) then 1586 call gettag(atlist,kk,natom,prtallatoms,tag) 1587 write(message,fmt)TRIM(message),ch10,thearray(:,kk),tag 1588 end if 1589 end do 1590 call wrtout(iout,message,'COLL') 1591 1592 end subroutine prtnatom
ABINIT/prtxfase [ Functions ]
NAME
prtxfase
FUNCTION
Print the values of xcart (X), forces (F) acell (A), Stresses (S), and energy (E) All values come from the history hist Also compute and print max and rms forces. Also compute absolute and relative differences with previous calculation
INPUTS
ab_mover<type abimover>=Subset of dtset only related with movement of ions and acell, contains: | dtion: Time step ! natom: Number of atoms | vis: viscosity | iatfix: Index of atoms and directions fixed | amass: Mass of ions hist<type abihist>=Historical record of positions, forces, stresses, cell and energies, itime= time step iout=unit number for printing
OUTPUT
(only writing)
SOURCE
1254 subroutine prtxfase(ab_mover,hist,itime,iout,pos) 1255 1256 !Arguments ------------------------------------ 1257 !scalars 1258 type(abimover),intent(in) :: ab_mover 1259 type(abihist),intent(in),target :: hist 1260 integer,intent(in) :: itime,iout 1261 integer,intent(in) :: pos 1262 !arrays 1263 1264 !Local variables------------------------------- 1265 !scalars 1266 integer :: jj,kk,unfixd,iprt 1267 real(dp) :: val_max,val_rms,ucvol ! Values maximal and RMS, Volume of Unitary cell 1268 real(dp) :: dEabs,dErel ! Diff of energy absolute and relative 1269 real(dp) :: ekin 1270 real(dp) :: angle(3),rmet(3,3) 1271 !character(len=80*(max(ab_mover%natom,3)+1)) :: msg 1272 !MGNAG: This is not very safe. One should use line-based output istead of appending chars 1273 ! and then outputting everything! For the time being I use this temporary hack to solve the problem with NAG 1274 character(len=max(80*(max(ab_mover%natom,3)+1),50000)) :: msg 1275 character(len=18) :: fmt1 1276 logical :: prtallatoms 1277 !arrays 1278 logical :: atlist(ab_mover%natom) 1279 real(dp),allocatable :: gred(:,:),xcart(:,:) 1280 real(dp),pointer :: acell(:),fcart(:,:),rprimd(:,:),strten(:),vel(:,:),xred(:,:) 1281 1282 ! *********************************************************** 1283 1284 fmt1='(a,a,1p,3e22.14)' 1285 1286 !########################################################## 1287 !### 1. Organize list of atoms to print 1288 1289 prtallatoms=.TRUE. 1290 do kk=1,ab_mover%natom 1291 if (ab_mover%prtatlist(kk)/=kk) prtallatoms=.FALSE. 1292 end do 1293 1294 atlist(:)=.FALSE. 1295 do iprt=1,ab_mover%natom 1296 if (ab_mover%prtatlist(iprt)>0.and.ab_mover%prtatlist(iprt)<=ab_mover%natom) atlist(ab_mover%prtatlist(iprt))=.TRUE. 1297 end do 1298 1299 acell => hist%acell(:,hist%ihist) 1300 rprimd => hist%rprimd(:,:,hist%ihist) 1301 xred => hist%xred(:,:,hist%ihist) 1302 fcart => hist%fcart(:,:,hist%ihist) 1303 strten => hist%strten(:,hist%ihist) 1304 vel => hist%vel(:,:,hist%ihist) 1305 1306 !########################################################### 1307 !### 1. Positions 1308 1309 ABI_MALLOC(xcart,(3,ab_mover%natom)) 1310 call xred2xcart(ab_mover%natom,rprimd,xcart,xred) 1311 1312 write(msg, '(a,a)' )ch10,' Cartesian coordinates (xcart) [bohr]' 1313 call prtnatom(atlist,iout,msg,ab_mover%natom,prtallatoms,xcart) 1314 1315 write(msg, '(a)' )' Reduced coordinates (xred)' 1316 call prtnatom(atlist,iout,msg,ab_mover%natom,prtallatoms,xred) 1317 1318 ABI_FREE(xcart) 1319 1320 !########################################################### 1321 !### 2. Forces 1322 1323 if(pos==mover_AFTER)then 1324 1325 ABI_MALLOC(gred,(3,ab_mover%natom)) 1326 call fcart2gred(fcart,gred,rprimd,ab_mover%natom) 1327 1328 ! Compute max |f| and rms f, 1329 ! EXCLUDING the components determined by iatfix 1330 val_max=0.0_dp 1331 val_rms=0.0_dp 1332 unfixd=0 1333 do kk=1,ab_mover%natom 1334 do jj=1,3 1335 if (ab_mover%iatfix(jj,kk) /= 1) then 1336 unfixd=unfixd+1 1337 val_rms=val_rms+fcart(jj,kk)**2 1338 val_max=max(val_max,abs(fcart(jj,kk)**2)) 1339 end if 1340 end do 1341 end do 1342 if ( unfixd /= 0 ) val_rms=sqrt(val_rms/dble(unfixd)) 1343 1344 write(msg, '(a,1p,2e12.5,a)' ) ' Cartesian forces (fcart) [Ha/bohr]; max,rms=',sqrt(val_max),val_rms,' (free atoms)' 1345 call prtnatom(atlist,iout,msg,ab_mover%natom,prtallatoms,fcart) 1346 1347 write(msg, '(a)' )' Gradient of E wrt nuclear positions in reduced coordinates (gred)' 1348 call prtnatom(atlist,iout,msg,ab_mover%natom,prtallatoms,gred) 1349 ABI_FREE(gred) 1350 end if 1351 1352 !########################################################### 1353 !### 3. Velocities 1354 1355 !Only if the velocities are being used 1356 if (hist%isVused)then 1357 ! Only if velocities are recorded in a history 1358 if (allocated(hist%vel))then 1359 ! Compute max |v| and rms v, 1360 ! EXCLUDING the components determined by iatfix 1361 val_max=0.0_dp 1362 val_rms=0.0_dp 1363 unfixd=0 1364 do kk=1,ab_mover%natom 1365 do jj=1,3 1366 if (ab_mover%iatfix(jj,kk) /= 1) then 1367 unfixd=unfixd+1 1368 val_rms=val_rms+vel(jj,kk)**2 1369 val_max=max(val_max,abs(vel(jj,kk)**2)) 1370 end if 1371 end do 1372 end do 1373 if ( unfixd /= 0 ) val_rms=sqrt(val_rms/dble(unfixd)) 1374 1375 write(msg, '(a,1p,2e12.5,a)' ) & 1376 & ' Cartesian velocities (vel) [bohr*Ha/hbar]; max,rms=',sqrt(val_max),val_rms,' (free atoms)' 1377 call prtnatom(atlist,iout,msg,ab_mover%natom,prtallatoms,vel) 1378 1379 ! Compute the ionic kinetic energy (no cell shape kinetic energy yet) 1380 ekin=0.0_dp 1381 do kk=1,ab_mover%natom 1382 do jj=1,3 1383 ! Warning : the fixing of atoms is implemented in reduced 1384 ! coordinates, so that this expression is wrong 1385 if (ab_mover%iatfix(jj,kk) == 0) then 1386 ekin=ekin+0.5_dp*ab_mover%amass(kk)*vel(jj,kk)**2 1387 end if 1388 end do 1389 end do 1390 write(msg, '(a,1p,e22.14,a)' )' Kinetic energy of ions (ekin) [Ha]=',ekin 1391 call wrtout(iout,msg,'COLL') 1392 end if 1393 end if 1394 1395 !########################################################### 1396 !### 3. ACELL 1397 1398 !Only if the acell is being used 1399 if (hist%isARused)then 1400 ! Only if acell is recorded in a history 1401 if (allocated(hist%acell))then 1402 write(msg, '(a)' )' Scale of Primitive Cell (acell) [bohr]' 1403 write(msg,fmt1)TRIM(msg),ch10,acell(:) 1404 call wrtout(iout,msg,'COLL') 1405 end if 1406 end if 1407 1408 !########################################################### 1409 !### 4. RPRIMD 1410 1411 !Only if the acell is being used 1412 if (hist%isARused)then 1413 ! Only if rprimd is recorded in a history 1414 if (allocated(hist%rprimd))then 1415 write(msg, '(a)' )' Real space primitive translations (rprimd) [bohr]' 1416 do kk=1,3 1417 write(msg,fmt1)TRIM(msg),ch10,rprimd(:,kk) 1418 end do 1419 call wrtout(iout,msg,'COLL') 1420 end if 1421 end if 1422 1423 !########################################################### 1424 !### 5. Unitary cell volume 1425 1426 if (ab_mover%optcell/=0)then 1427 1428 ucvol=& 1429 & rprimd(1,1)*(rprimd(2,2)*rprimd(3,3)-rprimd(3,2)*rprimd(2,3))+& 1430 & rprimd(2,1)*(rprimd(3,2)*rprimd(1,3)-rprimd(1,2)*rprimd(3,3))+& 1431 & rprimd(3,1)*(rprimd(1,2)*rprimd(2,3)-rprimd(2,2)*rprimd(1,3)) 1432 1433 write(msg, '(a,1p,e22.14)' )' Unitary Cell Volume (ucvol) [Bohr^3]=',ucvol 1434 call wrtout(iout,msg,'COLL') 1435 1436 ! ########################################################### 1437 ! ### 5. Angles and lengths 1438 1439 ! Compute real space metric. 1440 rmet = MATMUL(TRANSPOSE(rprimd),rprimd) 1441 1442 angle(1)=acos(rmet(2,3)/sqrt(rmet(2,2)*rmet(3,3)))/two_pi*360.0d0 1443 angle(2)=acos(rmet(1,3)/sqrt(rmet(1,1)*rmet(3,3)))/two_pi*360.0d0 1444 angle(3)=acos(rmet(1,2)/sqrt(rmet(1,1)*rmet(2,2)))/two_pi*360.0d0 1445 1446 write(msg, '(a)' )' Angles (23,13,12)= [degrees]' 1447 write(msg,fmt1)TRIM(msg),ch10,angle(:) 1448 call wrtout(iout,msg,'COLL') 1449 1450 write(msg, '(a)' ) ' Lengths [Bohr]' 1451 write(msg,fmt1)TRIM(msg),ch10,sqrt(rmet(1,1)),sqrt(rmet(2,2)),sqrt(rmet(3,3)) 1452 call wrtout(iout,msg,'COLL') 1453 1454 ! ########################################################### 1455 ! ### 5. Stress Tensor 1456 1457 if(pos==mover_AFTER)then 1458 ! Only if strten is recorded in a history 1459 if (allocated(hist%strten))then 1460 1461 write(msg, '(a)' ) ' Stress tensor in cartesian coordinates (strten) [Ha/bohr^3]' 1462 1463 write(msg,fmt1)TRIM(msg),ch10,strten(1),strten(6),strten(5) 1464 write(msg,fmt1)TRIM(msg),ch10,strten(6),strten(2),strten(4) 1465 write(msg,fmt1)TRIM(msg),ch10,strten(5),strten(4),strten(3) 1466 call wrtout(iout,msg,'COLL') 1467 end if 1468 end if 1469 end if 1470 1471 !########################################################### 1472 !### 6. Energy 1473 1474 if(pos==mover_AFTER)then 1475 write(msg, '(a,1p,e22.14)' )' Total energy (etotal) [Ha]=',hist%etot(hist%ihist) 1476 1477 if (itime>1)then 1478 jj = abihist_findIndex(hist,-1) 1479 dEabs=hist%etot(hist%ihist)-hist%etot(jj) 1480 dErel=2*dEabs/(abs(hist%etot(hist%ihist))+abs(hist%etot(jj))) 1481 write(msg, '(a,a,a,a)' )TRIM(msg),ch10,ch10,' Difference of energy with previous step (new-old):' 1482 write(msg, '(a,a,10a,a,1p,e12.5,a,10a,a,1p,e12.5)')& 1483 TRIM(msg),ch10,& 1484 (' ',jj=1,10),' Absolute (Ha)=',dEabs,ch10,& 1485 (' ',jj=1,10),' Relative =',dErel 1486 end if 1487 call wrtout(iout,msg,'COLL') 1488 end if 1489 1490 contains
ABINIT/wrt_moldyn_netcdf [ Functions ]
NAME
wrt_moldyn_netcdf
FUNCTION
Write two files for later molecular dynamics analysis: - MOLDYN.nc (netcdf format) : evolution of key quantities with time (pressure, energy, ...) - POSABIN : values of coordinates and velocities for the next time step
INPUTS
amass(natom)=mass of each atom, in unit of electronic mass (=amu*1822...) dtset <type(dataset_type)>=all input variables for this dataset itime=time step index option=1: write MOLDYN.nc file (netcdf format) 2: write POSABIN file 3: write both moldyn_file=name of the MD netcdf file mpi_enreg=information about MPI parallelization results_gs <type(results_gs_type)>=results (energy and its components, forces and its components, the stress tensor) of a ground-state computation rprimd(3,3)=real space primitive translations unpos=unit number for POSABIN file vel(3,natom)=velocities of atoms xred(3,natom)=reduced coordinates of atoms
OUTPUT
-- only printing --
SIDE EFFECTS
SOURCE
1630 subroutine wrt_moldyn_netcdf(amass,dtset,itime,option,moldyn_file,mpi_enreg,& 1631 & results_gs,rprimd,unpos,vel,xred) 1632 1633 use defs_basis 1634 use defs_abitypes 1635 use m_results_gs 1636 use m_abicore 1637 use m_errors 1638 #if defined HAVE_NETCDF 1639 use netcdf 1640 #endif 1641 1642 use m_io_tools, only : open_file, get_unit 1643 use m_geometry, only : xcart2xred, xred2xcart, metric 1644 1645 !Arguments ------------------------------------ 1646 !scalars 1647 integer,intent(in) :: itime,option,unpos 1648 character(fnlen),intent(in) :: moldyn_file 1649 type(dataset_type),intent(in) :: dtset 1650 type(MPI_type),intent(in) :: mpi_enreg 1651 type(results_gs_type),intent(in) :: results_gs 1652 !arrays 1653 real(dp),intent(in) :: amass(dtset%natom),rprimd(3,3) 1654 real(dp),intent(in),target :: vel(3,dtset%natom) 1655 real(dp),intent(in) :: xred(3,dtset%natom) 1656 1657 !Local variables------------------------------- 1658 !scalars 1659 integer,save :: ipos=0 1660 integer :: iatom,ii 1661 character(len=500) :: msg 1662 #if defined HAVE_NETCDF 1663 integer :: AtomNumDimid,AtomNumId,CelId,CellVolumeId,DimCoordid,DimScalarid,DimVectorid 1664 integer :: EkinDimid,EkinId,EpotDimid,EpotId,EntropyDimid,EntropyId,MassDimid,MassId,NbAtomsid 1665 integer :: ncerr,ncid,PosId,StressDimid,StressId,TensorSymDimid 1666 integer :: TimeDimid,TimestepDimid,TimestepId 1667 logical :: atom_fix 1668 real(dp) :: ekin,ucvol 1669 character(len=fnlen) :: ficname 1670 character(len=16) :: chain 1671 #endif 1672 !arrays 1673 #if defined HAVE_NETCDF 1674 integer :: PrimVectId(3) 1675 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3) 1676 real(dp),allocatable :: xcart(:,:) 1677 real(dp),pointer :: vcart(:,:),vred(:,:),vtmp(:,:) 1678 #endif 1679 1680 ! ************************************************************************* 1681 1682 !Only done by master processor, every nctime step 1683 if (mpi_enreg%me==0.and.dtset%nctime>0) then 1684 1685 ! Netcdf file name 1686 #if defined HAVE_NETCDF 1687 ficname = trim(moldyn_file)//'.nc' 1688 #endif 1689 1690 ! Xcart from Xred 1691 #if defined HAVE_NETCDF 1692 ABI_MALLOC(xcart,(3,dtset%natom)) 1693 call xred2xcart(dtset%natom,rprimd,xcart,xred) 1694 #endif 1695 1696 ! ========================================================================== 1697 ! First time step: write header of netcdf file 1698 ! ========================================================================== 1699 if (itime==1.and.(option==1.or.option==3)) then 1700 1701 ipos=0 1702 1703 #if defined HAVE_NETCDF 1704 ! Write message 1705 write(msg,'(4a)')ch10,' Open file ',trim(ficname),' to store molecular dynamics information.' 1706 call wrtout(std_out,msg,'COLL') 1707 1708 ! Create netcdf file 1709 ncerr = nf90_create(ficname, NF90_CLOBBER , ncid) 1710 NCF_CHECK_MSG(ncerr,'nf90_create') 1711 1712 ! Dimension time for netcdf (time dim is unlimited) 1713 ncerr = nf90_def_dim(ncid, "time", nf90_unlimited, TimeDimid) 1714 NCF_CHECK_MSG(ncerr,'nf90_def_dim') 1715 1716 ! Symetric Tensor Dimension 1717 ncerr = nf90_def_dim(ncid, "DimTensor", size(results_gs%strten), TensorSymDimid) 1718 NCF_CHECK_MSG(ncerr,'nf90_def_dim') 1719 1720 ! Coordinates Dimension 1721 ncerr = nf90_def_dim(ncid, "DimCoord", size(xcart,1), DimCoordid) 1722 NCF_CHECK_MSG(ncerr,'nf90_def_dim') 1723 1724 ! Atoms Dimensions 1725 ncerr = nf90_def_dim(ncid, "NbAtoms", dtset%natom, NbAtomsid) 1726 NCF_CHECK_MSG(ncerr,'nf90_def_dim') 1727 1728 ! Vector Dimension 1729 ncerr = nf90_def_dim(ncid, "DimVector", 3 , DimVectorid) 1730 NCF_CHECK_MSG(ncerr,'nf90_def_dim') 1731 1732 ! Scalar Dimension 1733 ncerr = nf90_def_dim(ncid, "DimScalar", 1 , DimScalarid) 1734 NCF_CHECK_MSG(ncerr,'nf90_def_dim') 1735 1736 ! Time step and time unit 1737 ncerr = nf90_def_var(ncid, "Time_step", nf90_double , DimScalarid, TimestepDimid) 1738 NCF_CHECK_MSG(ncerr,'nf90_def_var') 1739 ncerr = nf90_put_att(ncid, TimestepDimid, "units", "atomic time unit") 1740 NCF_CHECK_MSG(ncerr,'nf90_put_att') 1741 1742 ! Ionic masses 1743 ncerr = nf90_def_var(ncid, "Ionic_Mass", nf90_double , NbAtomsid, MassDimid) 1744 NCF_CHECK_MSG(ncerr,'nf90_def_var') 1745 ncerr = nf90_put_att(ncid, MassDimid, "units", "atomic mass unit") 1746 NCF_CHECK_MSG(ncerr,'nf90_put_att') 1747 1748 ! Ionic atomic numbers 1749 ncerr = nf90_def_var(ncid, "Ionic_Atomic_Number", nf90_double , NbAtomsid, AtomNumDimid) 1750 NCF_CHECK_MSG(ncerr,'nf90_def_var') 1751 1752 ! E_pot 1753 ncerr = nf90_def_var(ncid, "E_pot", nf90_double , TimeDimid, EpotDimid) 1754 NCF_CHECK_MSG(ncerr,'nf90_def_var') 1755 ncerr = nf90_put_att(ncid, EpotDimid, "units", "hartree") 1756 NCF_CHECK_MSG(ncerr,'nf90_put_att') 1757 1758 ! E_kin 1759 ncerr = nf90_def_var(ncid, "E_kin", nf90_double , TimeDimid, EkinDimid) 1760 NCF_CHECK_MSG(ncerr,'nf90_def_var') 1761 ncerr = nf90_put_att(ncid, EkinDimid, "units", "hartree") 1762 NCF_CHECK_MSG(ncerr,'nf90_put_att') 1763 1764 ! Entropy 1765 ncerr = nf90_def_var(ncid, "Entropy", nf90_double , TimeDimid, EntropyDimid) 1766 NCF_CHECK_MSG(ncerr,'nf90_def_var') 1767 ncerr = nf90_put_att(ncid, EntropyDimid, "units", "") 1768 NCF_CHECK_MSG(ncerr,'nf90_put_att') 1769 1770 ! Stress tensor 1771 ncerr = nf90_def_var(ncid, "Stress", nf90_double , (/TensorSymDimid,TimeDimid/), StressDimid) 1772 NCF_CHECK_MSG(ncerr,'nf90_def_var') 1773 ncerr = nf90_put_att(ncid, StressDimid, "units", "hartree/bohr^3") 1774 NCF_CHECK_MSG(ncerr,'nf90_put_att') 1775 1776 ! Positions 1777 ncerr = nf90_def_var(ncid, "Position", nf90_double ,(/DimCoordid,NbAtomsid,TimeDimid/), PosId) 1778 NCF_CHECK_MSG(ncerr,'nf90_def_var') 1779 ncerr = nf90_put_att(ncid, PosId, "units", "bohr") 1780 NCF_CHECK_MSG(ncerr,'nf90_put_att') 1781 1782 ! Celerities 1783 ncerr = nf90_def_var(ncid, "Celerity", nf90_double ,(/DimCoordid,NbAtomsid,TimeDimid/), CelId) 1784 NCF_CHECK_MSG(ncerr,'nf90_def_var') 1785 ncerr = nf90_put_att(ncid, CelId, "units", "bohr/(atomic time unit)") 1786 NCF_CHECK_MSG(ncerr,'nf90_put_att') 1787 1788 ! In case of volume cell constant 1789 if (dtset%optcell==0) then 1790 ! Primitive vectors 1791 do ii = 1,3 1792 write(unit=chain,fmt='(a15,i1)') "PrimitiveVector",ii 1793 ncerr = nf90_def_var(ncid, trim(chain), nf90_double , DimVectorid, PrimVectId(ii)) 1794 NCF_CHECK_MSG(ncerr,'nf90_def_var') 1795 end do 1796 ! Cell Volume 1797 ncerr = nf90_def_var(ncid, "Cell_Volume", nf90_double , DimScalarid, CellVolumeId) 1798 NCF_CHECK_MSG(ncerr,'nf90_def_var') 1799 ncerr = nf90_put_att(ncid, CellVolumeId, "units", "bohr^3") 1800 NCF_CHECK_MSG(ncerr,'nf90_put_att') 1801 end if 1802 1803 ! Leave define mode and close file 1804 ncerr = nf90_enddef(ncid) 1805 NCF_CHECK_MSG(ncerr,'nf90_enddef') 1806 ncerr = nf90_close(ncid) 1807 NCF_CHECK_MSG(ncerr,'nf90_close') 1808 #endif 1809 end if 1810 1811 ! ========================================================================== 1812 ! Write data to netcdf file (every nctime time step) 1813 ! ========================================================================== 1814 if (mod(itime, dtset%nctime)==0.and.(option==1.or.option==3)) then 1815 1816 ipos=ipos+1 1817 1818 #if defined HAVE_NETCDF 1819 ! Write message 1820 write(msg,'(3a)')ch10,' Store molecular dynamics information in file ',trim(ficname) 1821 call wrtout(std_out,msg,'COLL') 1822 1823 ! Open netcdf file 1824 ncerr = nf90_open(ficname, nf90_write, ncid) 1825 NCF_CHECK_MSG(ncerr,'nf90_open') 1826 1827 ! Time step 1828 ncerr = nf90_inq_varid(ncid, "Time_step", TimestepId) 1829 NCF_CHECK_MSG(ncerr,'nf90_inq_varid') 1830 ncerr = nf90_put_var(ncid, TimestepId, dtset%dtion) 1831 NCF_CHECK_MSG(ncerr,'nf90_put_var') 1832 1833 ! Ionic masses 1834 ncerr = nf90_inq_varid(ncid, "Ionic_Mass", MassId) 1835 NCF_CHECK_MSG(ncerr,'nf90_inq_varid') 1836 ncerr = nf90_put_var(ncid, MassId, amass, start = (/ 1 /), count=(/dtset%natom/)) 1837 NCF_CHECK_MSG(ncerr,'nf90_put_var') 1838 1839 ! Ionic atomic numbers 1840 ncerr = nf90_inq_varid(ncid, "Ionic_Atomic_Number", AtomNumId) 1841 NCF_CHECK_MSG(ncerr,'nf90_inq_varid') 1842 ncerr = nf90_put_var(ncid, AtomNumId, dtset%znucl(dtset%typat(:)),start=(/1/),count=(/dtset%natom/)) 1843 NCF_CHECK_MSG(ncerr,'nf90_put_var') 1844 1845 ! Epot 1846 ncerr = nf90_inq_varid(ncid, "E_pot", EpotId) 1847 NCF_CHECK_MSG(ncerr,'nf90_inq_varid') 1848 ncerr = nf90_put_var(ncid, EpotId, (/results_gs%etotal/), start=(/ipos/),count=(/1/)) 1849 NCF_CHECK_MSG(ncerr,'nf90_put_var') 1850 1851 ! Ekin 1852 ekin=zero;atom_fix=(maxval(dtset%iatfix)>0) 1853 if (dtset%ionmov==1.or.(.not.atom_fix)) then 1854 vcart => vel 1855 else 1856 ABI_MALLOC(vcart,(3,dtset%natom)) 1857 ABI_MALLOC(vred,(3,dtset%natom)) 1858 vtmp => vel 1859 call xcart2xred(dtset%natom,rprimd,vtmp,vred) 1860 do iatom=1,dtset%natom 1861 do ii=1,3 1862 if (dtset%iatfix(ii,iatom)==1) vred(ii,iatom)=zero 1863 end do 1864 end do 1865 call xred2xcart(dtset%natom,rprimd,vcart,vred) 1866 ABI_FREE(vred) 1867 end if 1868 do iatom=1,dtset%natom 1869 do ii=1,3 1870 ekin=ekin+half*amass(iatom)*vcart(ii,iatom)**2 1871 end do 1872 end do 1873 if (dtset%ionmov/=1.and.atom_fix) then 1874 ABI_FREE(vcart) 1875 end if 1876 ncerr = nf90_inq_varid(ncid, "E_kin", EkinId) 1877 NCF_CHECK_MSG(ncerr,'nf90_inq_varid') 1878 ncerr = nf90_put_var(ncid, EkinId, (/ekin/), start = (/ipos/),count=(/1/)) 1879 NCF_CHECK_MSG(ncerr,'nf90_put_var') 1880 1881 ! EntropyDimid 1882 ncerr = nf90_inq_varid(ncid, "Entropy", EntropyId) 1883 NCF_CHECK_MSG(ncerr,'nf90_inq_varid') 1884 ncerr = nf90_put_var(ncid, EntropyId, (/results_gs%energies%entropy/),start = (/ipos/),count=(/1/)) 1885 NCF_CHECK_MSG(ncerr,'nf90_put_var') 1886 1887 ! Stress tensor 1888 ncerr = nf90_inq_varid(ncid, "Stress", StressId) 1889 NCF_CHECK_MSG(ncerr,'nf90_inq_varid') 1890 ncerr = nf90_put_var(ncid, StressId, results_gs%strten, & 1891 & start=(/1,ipos/),count=(/size(results_gs%strten)/)) 1892 NCF_CHECK_MSG(ncerr,'nf90_put_var') 1893 1894 ! Positions 1895 ncerr = nf90_inq_varid(ncid, "Position", PosId) 1896 NCF_CHECK_MSG(ncerr,'nf90_inq_varid') 1897 ncerr = nf90_put_var(ncid, PosId, xcart, start=(/1,1,ipos/), & 1898 & count=(/size(xcart,1),dtset%natom,1/)) 1899 NCF_CHECK_MSG(ncerr,'nf90_put_var') 1900 1901 ! Celerities 1902 ncerr = nf90_inq_varid(ncid, "Celerity", CelId) 1903 NCF_CHECK_MSG(ncerr,'nf90_inq_varid') 1904 ncerr = nf90_put_var(ncid, CelId, vel, start=(/1,1,ipos/), & 1905 & count=(/size(vel,1),dtset%natom,1/) ) 1906 NCF_CHECK_MSG(ncerr,'nf90_put_var') 1907 1908 ! In case of volume cell constant 1909 if (dtset%optcell==0.and.ipos==1) then 1910 ! Primitive vectors 1911 do ii = 1,3 1912 write(unit=chain,fmt='(a15,i1)') "PrimitiveVector",ii 1913 ncerr = nf90_inq_varid(ncid, trim(chain), PrimVectId(ii) ) 1914 NCF_CHECK_MSG(ncerr,'nf90_inq_varid') 1915 ncerr = nf90_put_var(ncid, PrimVectId(ii), rprimd(:,ii)) 1916 NCF_CHECK_MSG(ncerr,'nf90_put_var') 1917 end do 1918 ! Cell Volume 1919 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 1920 ncerr = nf90_inq_varid(ncid, "Cell_Volume" , CellVolumeId) 1921 NCF_CHECK_MSG(ncerr,'nf90_inq_varid') 1922 ncerr = nf90_put_var(ncid, CellVolumeId, ucvol) 1923 NCF_CHECK_MSG(ncerr,'nf90_put_var') 1924 end if 1925 1926 ! Close file 1927 ncerr = nf90_close(ncid) 1928 NCF_CHECK_MSG(ncerr,'nf90_close') 1929 #endif 1930 end if 1931 1932 ! ========================================================================== 1933 ! Write data to POSABIN file (every nctime time step if option=3) 1934 ! ========================================================================== 1935 if ((mod(itime, dtset%nctime)==0.and.option==3).or.(option==2)) then 1936 1937 ! Open file for writing 1938 if (open_file('POSABIN',msg,unit=unpos,status='replace',form='formatted') /= 0 ) then 1939 ABI_ERROR(msg) 1940 end if 1941 1942 ! Write Positions 1943 if (dtset%natom>=1) write(unpos,'(a7,3d18.5)') 'xred ',(xred(ii,1),ii=1,3) 1944 if (dtset%natom>1) then 1945 do iatom=2,dtset%natom 1946 write(unpos,'(7x,3d18.5)') (xred(ii,iatom),ii=1,3) 1947 end do 1948 end if 1949 1950 ! Write Velocities 1951 if (dtset%natom>=1) write(unpos,'(a7,3d18.5)') 'vel ',(vel(ii,1),ii=1,3) 1952 if (dtset%natom>1) then 1953 do iatom=2,dtset%natom 1954 write(unpos,'(7x,3d18.5)') (vel(ii,iatom),ii=1,3) 1955 end do 1956 end if 1957 1958 ! Close file 1959 close(unpos) 1960 end if 1961 1962 #if defined HAVE_NETCDF 1963 ABI_FREE(xcart) 1964 #endif 1965 1966 ! ========================================================================== 1967 ! End if master proc 1968 end if 1969 1970 !Fake lines 1971 #if !defined HAVE_NETCDF 1972 if (.false.) write(std_out,*) moldyn_file,results_gs%etotal,rprimd(1,1) 1973 #endif 1974 1975 end subroutine wrt_moldyn_netcdf