TABLE OF CONTENTS
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
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (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 . For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
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=informations 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
PARENTS
mover
CHILDREN
metric,wrtout,xcart2xred,xred2xcart
SOURCE
47 #if defined HAVE_CONFIG_H 48 #include "config.h" 49 #endif 50 51 #include "abi_common.h" 52 53 subroutine wrt_moldyn_netcdf(amass,dtset,itime,option,moldyn_file,mpi_enreg,& 54 & results_gs,rprimd,unpos,vel,xred) 55 56 use defs_basis 57 use defs_abitypes 58 use m_results_gs 59 use m_profiling_abi 60 use m_errors 61 #if defined HAVE_NETCDF 62 use netcdf 63 #endif 64 65 use m_io_tools, only : open_file 66 use m_geometry, only : xcart2xred, xred2xcart, metric 67 68 !This section has been created automatically by the script Abilint (TD). 69 !Do not modify the following lines by hand. 70 #undef ABI_FUNC 71 #define ABI_FUNC 'wrt_moldyn_netcdf' 72 use interfaces_14_hidewrite 73 !End of the abilint section 74 75 implicit none 76 77 !Arguments ------------------------------------ 78 !scalars 79 integer,intent(in) :: itime,option,unpos 80 character(fnlen),intent(in) :: moldyn_file 81 type(dataset_type),intent(in) :: dtset 82 type(MPI_type),intent(in) :: mpi_enreg 83 type(results_gs_type),intent(in) :: results_gs 84 !arrays 85 real(dp),intent(in) :: amass(dtset%natom),rprimd(3,3) 86 real(dp),intent(in),target :: vel(3,dtset%natom) 87 real(dp),intent(in) :: xred(3,dtset%natom) 88 89 !Local variables------------------------------- 90 !scalars 91 integer,save :: ipos=0 92 integer :: iatom,ii 93 character(len=500) :: message 94 #if defined HAVE_NETCDF 95 integer :: AtomNumDimid,AtomNumId,CelId,CellVolumeId,DimCoordid,DimScalarid,DimVectorid 96 integer :: EkinDimid,EkinId,EpotDimid,EpotId,EntropyDimid,EntropyId,MassDimid,MassId,NbAtomsid 97 integer :: ncerr,ncid,PosId,StressDimid,StressId,TensorSymDimid 98 integer :: TimeDimid,TimestepDimid,TimestepId 99 logical :: atom_fix 100 real(dp) :: ekin,ucvol 101 character(len=fnlen) :: ficname 102 character(len=16) :: chain 103 #endif 104 !arrays 105 #if defined HAVE_NETCDF 106 integer :: PrimVectId(3) 107 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3) 108 real(dp),allocatable :: xcart(:,:) 109 real(dp),pointer :: vcart(:,:),vred(:,:),vtmp(:,:) 110 #endif 111 112 ! ************************************************************************* 113 114 !Only done by master processor, every nctime step 115 if (mpi_enreg%me==0.and.dtset%nctime>0) then 116 117 ! Netcdf file name 118 #if defined HAVE_NETCDF 119 ficname = trim(moldyn_file)//'.nc' 120 #endif 121 122 ! Xcart from Xred 123 #if defined HAVE_NETCDF 124 ABI_ALLOCATE(xcart,(3,dtset%natom)) 125 call xred2xcart(dtset%natom,rprimd,xcart,xred) 126 #endif 127 128 ! ========================================================================== 129 ! First time step: write header of netcdf file 130 ! ========================================================================== 131 if (itime==1.and.(option==1.or.option==3)) then 132 133 ipos=0 134 135 #if defined HAVE_NETCDF 136 ! Write message 137 write(message,'(4a)')ch10,' Open file ',trim(ficname),' to store molecular dynamics information.' 138 call wrtout(std_out,message,'COLL') 139 140 ! Create netcdf file 141 ncerr = nf90_create(ficname, NF90_CLOBBER , ncid) 142 NCF_CHECK_MSG(ncerr,'nf90_create') 143 144 ! Dimension time for netcdf (time dim is unlimited) 145 ncerr = nf90_def_dim(ncid, "time", nf90_unlimited, TimeDimid) 146 NCF_CHECK_MSG(ncerr,'nf90_def_dim') 147 148 ! Symetric Tensor Dimension 149 ncerr = nf90_def_dim(ncid, "DimTensor", size(results_gs%strten), TensorSymDimid) 150 NCF_CHECK_MSG(ncerr,'nf90_def_dim') 151 152 ! Coordinates Dimension 153 ncerr = nf90_def_dim(ncid, "DimCoord", size(xcart,1), DimCoordid) 154 NCF_CHECK_MSG(ncerr,'nf90_def_dim') 155 156 ! Atoms Dimensions 157 ncerr = nf90_def_dim(ncid, "NbAtoms", dtset%natom, NbAtomsid) 158 NCF_CHECK_MSG(ncerr,'nf90_def_dim') 159 160 ! Vector Dimension 161 ncerr = nf90_def_dim(ncid, "DimVector", 3 , DimVectorid) 162 NCF_CHECK_MSG(ncerr,'nf90_def_dim') 163 164 ! Scalar Dimension 165 ncerr = nf90_def_dim(ncid, "DimScalar", 1 , DimScalarid) 166 NCF_CHECK_MSG(ncerr,'nf90_def_dim') 167 168 ! Time step and time unit 169 ncerr = nf90_def_var(ncid, "Time_step", nf90_double , DimScalarid, TimestepDimid) 170 NCF_CHECK_MSG(ncerr,'nf90_def_var') 171 ncerr = nf90_put_att(ncid, TimestepDimid, "units", "atomic time unit") 172 NCF_CHECK_MSG(ncerr,'nf90_put_att') 173 174 ! Ionic masses 175 ncerr = nf90_def_var(ncid, "Ionic_Mass", nf90_double , NbAtomsid, MassDimid) 176 NCF_CHECK_MSG(ncerr,'nf90_def_var') 177 ncerr = nf90_put_att(ncid, MassDimid, "units", "atomic mass unit") 178 NCF_CHECK_MSG(ncerr,'nf90_put_att') 179 180 ! Ionic atomic numbers 181 ncerr = nf90_def_var(ncid, "Ionic_Atomic_Number", nf90_double , NbAtomsid, AtomNumDimid) 182 NCF_CHECK_MSG(ncerr,'nf90_def_var') 183 184 ! E_pot 185 ncerr = nf90_def_var(ncid, "E_pot", nf90_double , TimeDimid, EpotDimid) 186 NCF_CHECK_MSG(ncerr,'nf90_def_var') 187 ncerr = nf90_put_att(ncid, EpotDimid, "units", "hartree") 188 NCF_CHECK_MSG(ncerr,'nf90_put_att') 189 190 ! E_kin 191 ncerr = nf90_def_var(ncid, "E_kin", nf90_double , TimeDimid, EkinDimid) 192 NCF_CHECK_MSG(ncerr,'nf90_def_var') 193 ncerr = nf90_put_att(ncid, EkinDimid, "units", "hartree") 194 NCF_CHECK_MSG(ncerr,'nf90_put_att') 195 196 ! Entropy 197 ncerr = nf90_def_var(ncid, "Entropy", nf90_double , TimeDimid, EntropyDimid) 198 NCF_CHECK_MSG(ncerr,'nf90_def_var') 199 ncerr = nf90_put_att(ncid, EntropyDimid, "units", "") 200 NCF_CHECK_MSG(ncerr,'nf90_put_att') 201 202 ! Stress tensor 203 ncerr = nf90_def_var(ncid, "Stress", nf90_double , (/TensorSymDimid,TimeDimid/), StressDimid) 204 NCF_CHECK_MSG(ncerr,'nf90_def_var') 205 ncerr = nf90_put_att(ncid, StressDimid, "units", "hartree/bohr^3") 206 NCF_CHECK_MSG(ncerr,'nf90_put_att') 207 208 ! Positions 209 ncerr = nf90_def_var(ncid, "Position", nf90_double ,(/DimCoordid,NbAtomsid,TimeDimid/), PosId) 210 NCF_CHECK_MSG(ncerr,'nf90_def_var') 211 ncerr = nf90_put_att(ncid, PosId, "units", "bohr") 212 NCF_CHECK_MSG(ncerr,'nf90_put_att') 213 214 ! Celerities 215 ncerr = nf90_def_var(ncid, "Celerity", nf90_double ,(/DimCoordid,NbAtomsid,TimeDimid/), CelId) 216 NCF_CHECK_MSG(ncerr,'nf90_def_var') 217 ncerr = nf90_put_att(ncid, CelId, "units", "bohr/(atomic time unit)") 218 NCF_CHECK_MSG(ncerr,'nf90_put_att') 219 220 ! In case of volume cell constant 221 if (dtset%optcell==0) then 222 ! Primitive vectors 223 do ii = 1,3 224 write(unit=chain,fmt='(a15,i1)') "PrimitiveVector",ii 225 ncerr = nf90_def_var(ncid, trim(chain), nf90_double , DimVectorid, PrimVectId(ii)) 226 NCF_CHECK_MSG(ncerr,'nf90_def_var') 227 end do 228 ! Cell Volume 229 ncerr = nf90_def_var(ncid, "Cell_Volume", nf90_double , DimScalarid, CellVolumeId) 230 NCF_CHECK_MSG(ncerr,'nf90_def_var') 231 ncerr = nf90_put_att(ncid, CellVolumeId, "units", "bohr^3") 232 NCF_CHECK_MSG(ncerr,'nf90_put_att') 233 end if 234 235 ! Leave define mode and close file 236 ncerr = nf90_enddef(ncid) 237 NCF_CHECK_MSG(ncerr,'nf90_enddef') 238 ncerr = nf90_close(ncid) 239 NCF_CHECK_MSG(ncerr,'nf90_close') 240 #endif 241 end if 242 243 ! ========================================================================== 244 ! Write data to netcdf file (every nctime time step) 245 ! ========================================================================== 246 if (mod(itime, dtset%nctime)==0.and.(option==1.or.option==3)) then 247 248 ipos=ipos+1 249 250 #if defined HAVE_NETCDF 251 ! Write message 252 write(message,'(3a)')ch10,' Store molecular dynamics information in file ',trim(ficname) 253 call wrtout(std_out,message,'COLL') 254 255 ! Open netcdf file 256 ncerr = nf90_open(ficname, nf90_write, ncid) 257 NCF_CHECK_MSG(ncerr,'nf90_open') 258 259 ! Time step 260 ncerr = nf90_inq_varid(ncid, "Time_step", TimestepId) 261 NCF_CHECK_MSG(ncerr,'nf90_inq_varid') 262 ncerr = nf90_put_var(ncid, TimestepId, dtset%dtion) 263 NCF_CHECK_MSG(ncerr,'nf90_put_var') 264 265 ! Ionic masses 266 ncerr = nf90_inq_varid(ncid, "Ionic_Mass", MassId) 267 NCF_CHECK_MSG(ncerr,'nf90_inq_varid') 268 ncerr = nf90_put_var(ncid, MassId, amass, start = (/ 1 /), count=(/dtset%natom/)) 269 NCF_CHECK_MSG(ncerr,'nf90_put_var') 270 271 ! Ionic atomic numbers 272 ncerr = nf90_inq_varid(ncid, "Ionic_Atomic_Number", AtomNumId) 273 NCF_CHECK_MSG(ncerr,'nf90_inq_varid') 274 ncerr = nf90_put_var(ncid, AtomNumId, dtset%znucl(dtset%typat(:)),start=(/1/),count=(/dtset%natom/)) 275 NCF_CHECK_MSG(ncerr,'nf90_put_var') 276 277 ! Epot 278 ncerr = nf90_inq_varid(ncid, "E_pot", EpotId) 279 NCF_CHECK_MSG(ncerr,'nf90_inq_varid') 280 ncerr = nf90_put_var(ncid, EpotId, (/results_gs%etotal/), start=(/ipos/),count=(/1/)) 281 NCF_CHECK_MSG(ncerr,'nf90_put_var') 282 283 ! Ekin 284 ekin=zero;atom_fix=(maxval(dtset%iatfix)>0) 285 if (dtset%ionmov==1.or.(.not.atom_fix)) then 286 vcart => vel 287 else 288 ABI_ALLOCATE(vcart,(3,dtset%natom)) 289 ABI_ALLOCATE(vred,(3,dtset%natom)) 290 vtmp => vel 291 call xcart2xred(dtset%natom,rprimd,vtmp,vred) 292 do iatom=1,dtset%natom 293 do ii=1,3 294 if (dtset%iatfix(ii,iatom)==1) vred(ii,iatom)=zero 295 end do 296 end do 297 call xred2xcart(dtset%natom,rprimd,vcart,vred) 298 ABI_DEALLOCATE(vred) 299 end if 300 do iatom=1,dtset%natom 301 do ii=1,3 302 ekin=ekin+half*amass(iatom)*vcart(ii,iatom)**2 303 end do 304 end do 305 if (dtset%ionmov/=1.and.atom_fix) then 306 ABI_DEALLOCATE(vcart) 307 end if 308 ncerr = nf90_inq_varid(ncid, "E_kin", EkinId) 309 NCF_CHECK_MSG(ncerr,'nf90_inq_varid') 310 ncerr = nf90_put_var(ncid, EkinId, (/ekin/), start = (/ipos/),count=(/1/)) 311 NCF_CHECK_MSG(ncerr,'nf90_put_var') 312 313 ! EntropyDimid 314 ncerr = nf90_inq_varid(ncid, "Entropy", EntropyId) 315 NCF_CHECK_MSG(ncerr,'nf90_inq_varid') 316 ncerr = nf90_put_var(ncid, EntropyId, (/results_gs%energies%entropy/),start = (/ipos/),count=(/1/)) 317 NCF_CHECK_MSG(ncerr,'nf90_put_var') 318 319 ! Stress tensor 320 ncerr = nf90_inq_varid(ncid, "Stress", StressId) 321 NCF_CHECK_MSG(ncerr,'nf90_inq_varid') 322 ncerr = nf90_put_var(ncid, StressId, results_gs%strten, & 323 & start=(/1,ipos/),count=(/size(results_gs%strten)/)) 324 NCF_CHECK_MSG(ncerr,'nf90_put_var') 325 326 ! Positions 327 ncerr = nf90_inq_varid(ncid, "Position", PosId) 328 NCF_CHECK_MSG(ncerr,'nf90_inq_varid') 329 ncerr = nf90_put_var(ncid, PosId, xcart, start=(/1,1,ipos/), & 330 & count=(/size(xcart,1),dtset%natom,1/)) 331 NCF_CHECK_MSG(ncerr,'nf90_put_var') 332 333 ! Celerities 334 ncerr = nf90_inq_varid(ncid, "Celerity", CelId) 335 NCF_CHECK_MSG(ncerr,'nf90_inq_varid') 336 ncerr = nf90_put_var(ncid, CelId, vel, start=(/1,1,ipos/), & 337 & count=(/size(vel,1),dtset%natom,1/) ) 338 NCF_CHECK_MSG(ncerr,'nf90_put_var') 339 340 ! In case of volume cell constant 341 if (dtset%optcell==0.and.ipos==1) then 342 ! Primitive vectors 343 do ii = 1,3 344 write(unit=chain,fmt='(a15,i1)') "PrimitiveVector",ii 345 ncerr = nf90_inq_varid(ncid, trim(chain), PrimVectId(ii) ) 346 NCF_CHECK_MSG(ncerr,'nf90_inq_varid') 347 ncerr = nf90_put_var(ncid, PrimVectId(ii), rprimd(:,ii)) 348 NCF_CHECK_MSG(ncerr,'nf90_put_var') 349 end do 350 ! Cell Volume 351 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 352 ncerr = nf90_inq_varid(ncid, "Cell_Volume" , CellVolumeId) 353 NCF_CHECK_MSG(ncerr,'nf90_inq_varid') 354 ncerr = nf90_put_var(ncid, CellVolumeId, ucvol) 355 NCF_CHECK_MSG(ncerr,'nf90_put_var') 356 end if 357 358 ! Close file 359 ncerr = nf90_close(ncid) 360 NCF_CHECK_MSG(ncerr,'nf90_close') 361 #endif 362 end if 363 364 ! ========================================================================== 365 ! Write data to POSABIN file (every nctime time step if option=3) 366 ! ========================================================================== 367 if ((mod(itime, dtset%nctime)==0.and.option==3).or.(option==2)) then 368 369 ! Open file for writing 370 if (open_file('POSABIN',message,unit=unpos,status='replace',form='formatted') /= 0 ) then 371 MSG_ERROR(message) 372 end if 373 374 ! Write Positions 375 if (dtset%natom>=1) write(unpos,'(a7,3d18.5)') 'xred ',(xred(ii,1),ii=1,3) 376 if (dtset%natom>1) then 377 do iatom=2,dtset%natom 378 write(unpos,'(7x,3d18.5)') (xred(ii,iatom),ii=1,3) 379 end do 380 end if 381 382 ! Write Velocities 383 if (dtset%natom>=1) write(unpos,'(a7,3d18.5)') 'vel ',(vel(ii,1),ii=1,3) 384 if (dtset%natom>1) then 385 do iatom=2,dtset%natom 386 write(unpos,'(7x,3d18.5)') (vel(ii,iatom),ii=1,3) 387 end do 388 end if 389 390 ! Close file 391 close(unpos) 392 end if 393 394 #if defined HAVE_NETCDF 395 ABI_DEALLOCATE(xcart) 396 #endif 397 398 ! ========================================================================== 399 ! End if master proc 400 end if 401 402 !Fake lines 403 #if !defined HAVE_NETCDF 404 if (.false.) write(std_out,*) moldyn_file,results_gs%etotal,rprimd(1,1) 405 #endif 406 407 end subroutine wrt_moldyn_netcdf