forum@abinit.org
Subject: The ABINIT Users Mailing List ( CLOSED )
List archive
- From: Matthieu Verstraete <mjv500@york.ac.uk>
- To: forum@abinit.org
- Subject: Re: [abinit-forum] Problem with cut3d and Xcrysden format
- Date: Fri, 6 Jul 2007 15:51:38 +0100 (BST)
Hello Mauro and co.
the fortran format was insufficient for your numbers. I'm surprised this hasn't happened before. Anyhow I have modified the source and committed to 5.4.3, which should appear soon.
For the impatient, use the attached cut3d.F90
Matthieu
--
================================================================
Dr. Matthieu Verstraete mailto:mjv500@york.ac.uk
Dept. of Physics, University of York, tel: +44 1904 43 22 08
Heslington, YO10 5DD York, United Kingdom fax: +44 1904 43 22 14!{\src2tex{textfont=tt}}
!!****p* ABINIT/cut3d
!! NAME
!! cut3d
!!
!! FUNCTION
!! Main routine for the analysis of the density and potential files,
!! as well as other files with the ABINIT header.
!!
!! COPYRIGHT
!! Copyright (C) 1999-2007 ABINIT group (GMR, RC, LSI, XG, NCJ, JFB, MCote,
LPizzagalli)
!! 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
!! (main program)
!!
!! OUTPUT
!! (main program)
!!
!! PARENTS
!!
!! CHILDREN
!! date_and_time,handle_ncerr,hdr_clean,hdr_io,herald,hirsh,lineint
!! localorb_s,planeint,pointint,rrho,rtau,volumeint,wffile
!!
!! SOURCE
#if defined HAVE_CONFIG_H
#include "config.h"
#endif
program cut3d
use defs_basis
#if defined FC_NAG
use f90_unix
#endif
#if defined HAVE_NETCDF
use netcdf
#endif
use defs_datatypes
use defs_infos
!This section has been created automatically by the script Abilint (TD). Do
not modify these by hand.
#ifdef HAVE_FORTRAN_INTERFACES
use interfaces_01manage_mpi
use interfaces_13io_mpi
use interfaces_14iowfdenpot
use interfaces_19cut3d
#else
use defs_interfaces
#endif
!End of the abilint section
implicit none
!Arguments -----------------------------------
!Local variables-------------------------------
! natom = number of atoms in the unit cell
! nr1,nr2,nr3 = grid size (nr1 x nr2 x nr3 = filrho dimension)
! ntypat = number of atom types
! ucvol = unit cell volume (> 0)
! densfileformat = flag for the format of the density file:
! 0 = ASCII
! 1 = binary
! denval = density value exported by interpol, to be wrote in the output
file
! filrho = name of the density file (ASCII or binary)
! filtau = name of the atomic position file (Xmol format)
! acell = unit cell parameters
! rprim = orientation of the unit cell axes
character(len=*), parameter :: INPUTfile='cut.in'
!For NetCDF *******************************************************
#if defined HAVE_NETCDF
character(len=3), parameter ::
monnam(12)=(/'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'/)
#endif
!scalars
integer :: atomicnumVarID,atomposiVarID,dataVarID,dd,densfileformat,exchn2n3
integer :: fform0,grid1VarID,grid2VarID,grid3VarID,gridshift1,gridshift2
integer :: gridshift3,gridsize1DimID,gridsize2DimID,gridsize3DimID,i1,i2,i3
integer :: iatom,igrid,ii,ii1,ii2,ii3,index,iprompt,ir1,ir2,ir3,ispden,itask
integer :: latDimID,latticevecVarID,mm,natom,nbatomDimID,ncid,nr1,nr2,nr3
integer :: nspden,ntypat,originVarID,posDimID,rdwr,status,titlechoice,unitfi
integer :: yyyy
real(dp) :: denval,ucvol
real*4 :: xm,xnow,xp,ym,ynow,yp,zm,znow,zp
logical :: filexist
character(1) :: outputchar
character(len=10) :: strtime
character(len=11) :: stridate
character(len=24) :: codename
character(len=5) :: strzone
character(len=50) :: chr_inputfname
character(len=500) :: filetitle
character(len=8) :: strdat
character(len=fnlen) :: filnam,filrho,filtau
type(hdr_type) :: hdr
!arrays
integer :: values(8)
real(dp) :: acell(3),gridwavefun1(3,2),gridwavefun2(3,2),gridwavefun3(3,2)
real(dp) :: originatt(3,3),rprim(3,3),rprimd(3,3),shift_tau(3),value(2)
real(dp) :: xcart2(3)
real(dp),allocatable :: grid(:,:,:),grid_full(:,:,:,:),nuclz(:),tau2(:,:)
!DEBUG (LD)
real(dp),allocatable ::
gridtt(:,:,:),gridux(:,:,:),griddy(:,:,:),gridmz(:,:,:)
!ENDDEBUG (LD)
real(dp),allocatable :: xcart(:,:),xred(:,:)
real*4,allocatable :: rhomacu(:,:)
!******************************************************************
!BEGIN EXECUTABLE SECTION
codename='CUT3D '//repeat(' ',18)
call herald(codename,abinit_version,std_out)
!Get name of density file
write(*,*)
write(*,*) ' What is the name of the 3D function (density, potential or
wavef) file ?'
read(5,*)filrho
write(*,*) ' => Your 3D function file is : ',trim(filrho)
write(*,*)
!Checking the existence of data file
inquire (file=filrho,exist=filexist)
if (.NOT. filexist) then
write(*,*) 'Error, missing data file: ',filrho
stop
end if
!Get its type
write(*,*) ' Does this file contain formatted 3D ASCII data (=0) '
write(*,*) ' or unformatted binary header + 3D data (=1) ?'
read(5,*)densfileformat
!Treat the different cases : formatted or unformatted
if(densfileformat==1)then
write(*,*) ' 1 => Your file contains unformatted binary header + 3D data'
write(*,*) ' The information it contains should be sufficient.'
open(unit=19,file=filrho,form='unformatted',status='old')
write(6, '(a,a,a,i4)' )&
& ' cut3d : read file ',trim(filrho),' from unit number 19.'
write(*,*)
! Read the header
rdwr=1 ; unitfi=19
call hdr_io(fform0,hdr,rdwr,unitfi)
! Echo part of the header
rdwr=4 ; unitfi=6
call hdr_io(fform0,hdr,rdwr,unitfi)
natom=hdr%natom
nr1=hdr%ngfft(1)
nr2=hdr%ngfft(2)
nr3=hdr%ngfft(3)
nspden=hdr%nspden
ntypat=hdr%ntypat
rprimd(:,:)=hdr%rprimd(:,:)
! Need to know natom in order to allocate xcart
allocate(xcart(3,natom),xred(3,natom))
xred(:,:)=hdr%xred(:,:)
do iatom=1,natom
xcart(:,iatom)=rprimd(:,1)*xred(1,iatom)+&
& rprimd(:,2)*xred(2,iatom)+&
& rprimd(:,3)*xred(3,iatom)
end do
if(fform0>50)then
ispden=0
if(nspden/=1)then
write(*, '(a)' )' '
write(*, '(a)' )' * This file contains more than one spin component,'
write(*, '(a,i3,a)' )' (indeed, nspden=',nspden,' )'
write(*, '(a)' )' Some of the tasks that you will define later will
concern all spin components.'
write(*, '(a)' )' Others tasks might require you to have chosen among
the following :'
endif
if(nspden==2)then
write(*, '(a)' )' ispden= 0 ==> Total density'
write(*, '(a)' )' ispden= 1 ==> spin-up density'
write(*, '(a)' )' ispden= 2 ==> spin-down density'
write(*, '(a)' )' ispden= 3 ==> spin-polarization (or magnetization)
density'
write(*, '(a)' )' spin up - spin down difference.'
endif
if(nspden==4)then
write(*, '(a)' )' ispden= 0 ==> Total density'
write(*, '(a)' )' ispden= 1 ==> magnetization in the x direction'
write(*, '(a)' )' ispden= 2 ==> magnetization in the y direction'
write(*, '(a)' )' ispden= 3 ==> magnetization in the z direction'
write(*, '(a)' )' ispden= 4 might be used to plot the magnetization
(3D) in the XCrysDen format,'
endif
if(nspden/=1)then
write(*,*)' Please define ispden :'
read(5,*)ispden
write(6, '(a,i3)' )' You entered ispden=',ispden
end if
end if
else if(densfileformat==0)then
write(*,*) ' 0 => Your file contains formatted 3D ASCII data'
write(*,*) ' The complementary information is taken from ',trim(INPUTfile)
! Checking the existence of input file
inquire (file=INPUTfile,exist=filexist)
if (.NOT. filexist) then
write(*,*) 'Error, missing input file ',INPUTfile
stop
end if
! Read in the input file INPUTfile
write(*,*)
write(*,*) 'READING FROM FILE ',INPUTfile
open(32,file=INPUTfile,status='old')
read(32,'(a)') filtau
! Checking the existence of atomic position file
inquire (file=filtau,exist=filexist)
if (.NOT. filexist) then
write(*,*) 'Error, missing atomic position file ',trim(filtau)
stop
end if
write(*,*) 'atomic position file (Xmol format):',trim(filtau)
! Read cell parameters
read(32,*) acell(1),acell(2),acell(3)
do ii=1,3
if (acell(ii) <= 0.0) then
write(*,*) 'Error, invalid value for acell(',ii,'): ',acell(ii)
stop
end if
end do
read(32,*) rprim
do ii=1,3
rprimd(:,ii)=rprim(:,ii)*acell(ii)
end do
! FFT grid, number of atoms, number of atom types
read(32,*) nr1,nr2,nr3
read(32,*) natom
read(32,*) ntypat
close(32)
! Need to know natom in order to allocate xcart
allocate(xcart(3,natom))
call rtau(filtau,xcart,natom,ntypat)
write(*,*)
! By default there is only one spin component, and one works with a total
density
nspden=1 ; ispden=0 ; fform0=52
open(unit=19,file=filrho,form='formatted',status='old')
else
write(*,*)'Error, value for density file format is invalid: ',densfileformat
stop
end if
write(*,*)
write(*,*) '==========================================================='
write(*,*)
!Echo the value of different input parameters
write(*,*)'ECHO important input variables ...'
write(*,*)
write(*,*) ' Dimensional primitive vectors (ABINIT equivalent : rprimd):'
write(*, '(3es16.6)' ) rprimd(1:3,1)
write(*, '(3es16.6)' ) rprimd(1:3,2)
write(*, '(3es16.6)' ) rprimd(1:3,3)
!Here test the non-colinearity of rprim vectors
ucvol=rprimd(1,1)*(rprimd(2,2)*rprimd(3,3)-rprimd(3,2)*rprimd(2,3))+&
& rprimd(2,1)*(rprimd(3,2)*rprimd(1,3)-rprimd(1,2)*rprimd(3,3))+&
& rprimd(3,1)*(rprimd(1,2)*rprimd(2,3)-rprimd(2,2)*rprimd(1,3))
if (abs(ucvol)<1.0d-12) then
write(*,*)' At least two rprimd(,) vectors are colinear'
write(*,*)' Please check the input rprim and acell, or rprimd.'
write(*,*)' The program will stop'
stop
end if
write(*, '(a,3i5)' ) ' Grid density (ABINIT equivalent : ngfft):
',nr1,nr2,nr3
write(*,*) ' Number of atoms :',natom
write(*,*) ' Number of atomic types:',ntypat
write(*,*)
write(*,*) ' # Atomic positions (cartesian coordinates - Bohr)'
do iatom=1,natom
write(6, '(i4,3es16.6)' )iatom,xcart(1:3,iatom)
end do
write(*,*)
!------------------------------------------------------------------------
!Branching : either WF file, or DEN/POT file.
if(densfileformat==1 .and. fform0<50)then
write(*,*)' This file is a WF file. '
exchn2n3=0
write(*,*)" If you want to analyze one wavefunction, type
0 "
write(*,*)" If you want to construct Wannier-type Localized Orbitals, type
2 "
read(*,*)ii1
write(*,*)" You typed ",ii1
if(ii1==0)then
call wffile(hdr%ecut_eff,exchn2n3,hdr%headform,hdr%istwfk,hdr%kptns,&
& natom,hdr%nband,hdr%nkpt,hdr%npwarr,&
&
nr1,nr2,nr3,hdr%nspinor,hdr%nsppol,ntypat,rprimd,xcart,hdr%typat,hdr%znucltypat)
else if(ii1==2)then
! Read the name of the input file name :
read(5,'(a)')chr_inputfname
call
localorb_S(chr_inputfname,hdr%ecut_eff,exchn2n3,hdr%headform,hdr%istwfk,hdr%kptns,&
& natom,hdr%nband,hdr%nkpt,hdr%npwarr,&
&
nr1,nr2,nr3,hdr%nspinor,hdr%nsppol,ntypat,rprimd,xcart,hdr%typat,hdr%znucltypat)
write(*,*)" "
write(*,*)" ###################################################### "
write(*,*)" "
write(*,*)" Localized orbital files fort.1**1 for spin up "
write(*,*)" and fort.1**2 for spin dn written."
write(*,*)" "
write(*,*)" ###################################################### "
else
write(*,*)" Option ",ii1," is not allowed => stop "
end if
!-------------------------------------------------------------------------
else ! This is a DEN/POT file
! This should become a subroutine
write(*,*)' This file is a Density or Potential file '
! Read the function on the 3D grid
allocate(grid(nr1,nr2,nr3),grid_full(nr1,nr2,nr3,nspden))
allocate(gridtt(nr1,nr2,nr3),gridux(nr1,nr2,nr3))
allocate(griddy(nr1,nr2,nr3),gridmz(nr1,nr2,nr3))
call rrho(densfileformat,grid_full,nr1,nr2,nr3,nspden)
! Do not forget that the first sub-array of a density file is the total
density,
! while the first sub-array of a potential file is the spin-up potential
if(fform0==51 .or. fform0==52)then ! Density case
! gridtt= grid --> Total density.
! gridux= grid --> spin-Up density, or magnetization density in X direction.
! griddy= grid --> spin-Down density, or magnetization density in Y
direction.
! gridmz= grid --> spin-polarization density (Magnetization),
! or magnetization density in Z direction.
gridtt(:,:,:)=grid_full(:,:,:,1)
if(nspden==2)then
gridux(:,:,:)=grid_full(:,:,:,2)
griddy(:,:,:)=grid_full(:,:,:,1)-grid_full(:,:,:,2)
gridmz(:,:,:)=-grid_full(:,:,:,1)+two*grid_full(:,:,:,2)
else if(nspden==4)then
gridux(:,:,:)=grid_full(:,:,:,2)
griddy(:,:,:)=grid_full(:,:,:,3)
gridmz(:,:,:)=grid_full(:,:,:,4)
endif
if(nspden==1)then
grid(:,:,:)=grid_full(:,:,:,1)
else
if(ispden==0)then
grid(:,:,:)=gridtt(:,:,:)
else if(ispden==1)then
grid(:,:,:)=gridux(:,:,:)
else if(ispden==2)then
grid(:,:,:)=griddy(:,:,:)
else if(ispden==3)then
grid(:,:,:)=gridmz(:,:,:)
! if(ispden==0)then
! grid(:,:,:)=grid_full(:,:,:,1)
! else if(ispden==1)then
! grid(:,:,:)=grid_full(:,:,:,2)
! else if(ispden==2)then
! grid(:,:,:)=grid_full(:,:,:,1)-grid_full(:,:,:,2)
! else if(ispden==-1)then
! grid(:,:,:)=-grid_full(:,:,:,1)+two*grid_full(:,:,:,2)
else if(ispden==4)then
write(*,*) ' '
else
write(*,*)' Error: bad ispden value'
stop
end if
endif
else if(fform0==101 .or. fform0==102)then ! Potential case
if(ispden==0)then
grid(:,:,:)=grid_full(:,:,:,1)
else if(ispden==1 .or. ispden==2)then
grid(:,:,:)=grid_full(:,:,:,ispden)
else
write(*,*)' Error: bad ispden value'
stop
end if
end if
write(*,*)
write(*,*) ' 3D function was read. Ready for further treatment.'
write(*,*)
write(*,*) '==========================================================='
write(*,*)
!------------------------------------------------------------------------
! At this moment all the input is done
! The code knows the geometry of the system,
! and the data file (electron density, potential, etc).
! It will further calculate the electron density by interpolation in
! a point, along a line or in a plane.
do
do
write(*,*) ' What is your choice ? Type:'
write(*,*) ' 0 => exit'
write(*,*) ' 1 => point (interpolation of data for a single point)'
write(*,*) ' 2 => line (interpolation of data along a line)'
write(*,*) ' 3 => plane (interpolation of data in a plane)'
write(*,*) ' 4 => volume (interpolation of data in a volume)'
write(*,*) ' 5 => 3D formatted data (output the bare 3D data - one
column)'
write(*,*) ' 6 => 3D indexed data (bare 3D data, preceeded by 3D index)'
write(*,*) ' 7 => 3D Molekel formatted data '
write(*,*) ' 8 => 3D data with coordinates (tecplot ASCII format)'
write(*,*) ' 9 => output .xsf file for XCrysDen'
write(*,*) ' 10 => output .dx file for OpenDx'
write(*,*) ' 11 => compute atomic charge using the Hirshfeld method'
write(*,*) ' 12 => NetCDF file'
read(*,*) itask
write(*, '(a,a,i2,a)' ) ch10,' Your choice is ',itask,ch10
if( 5<=itask .and. itask<=10 .or. itask==12 )then
write(*,*) ch10,' Enter the name of an output file:'
read(*,*) filnam
write(*,*) ' The name of your file is : ',trim(filnam)
end if
select case(itask)
!DEBUG (LD)
case(1) ! point calculation
! call pointint(grid,nr1,nr2,nr3,rprimd)
call
pointint(grid,gridtt,gridux,griddy,gridmz,nr1,nr2,nr3,nspden,rprimd)
exit
case(2) ! line calculation
! call lineint(grid,nr1,nr2,nr3,rprimd)
call lineint(grid,gridtt,gridux,griddy,gridmz,nr1,nr2,nr3,nspden,rprimd)
exit
case(3) ! plane calculation
! call planeint(grid,natom,nr1,nr2,nr3,rprimd,xcart)
call
planeint(grid,gridtt,gridux,griddy,gridmz,natom,nr1,nr2,nr3,nspden,rprimd,xcart)
exit
case(4) ! volume calculation
write(*,*) ' Enter volume calculation'
! call volumeint(grid,natom,nr1,nr2,nr3,rprimd,xcart)
call
volumeint(grid,gridtt,gridux,griddy,gridmz,natom,nr1,nr2,nr3,nspden,rprimd,xcart)
exit
case(5)
! Rewrite the data on a formatted file, just in one (or four) column(s)
open(unit=31,file=trim(filnam),status='unknown')
if(nspden==1)then
do i3=1,nr3
do i2=1,nr2
do i1=1,nr1
write(31,'(4(es22.12))') grid(i1,i2,i3)
end do
end do
end do
else
do i3=1,nr3
do i2=1,nr2
do i1=1,nr1
write(31,'(4(es22.12))') gridtt(i1,i2,i3), gridux(i1,i2,i3),
griddy(i1,i2,i3), gridmz(i1,i2,i3)
end do
end do
end do
endif
close(31)
exit
case(6)
! Rewrite the data on a formatted file, 3D index + density
open(unit=31,file=trim(filnam),status='unknown')
if(nspden==1)then
write(31,*)' i1 i2 i3 data '
do i3=1,nr3
do i2=1,nr2
do i1=1,nr1
write(31,'(3i6,4(es24.14))') i1,i2,i3,grid(i1,i2,i3)
end do
end do
end do
else
if(nspden==2)then
write(31,*)' i1 i2 i3 non-spin-polarized spin up spin
down difference '
else if(nspden==4)then
write(31,*)' i1 i2 i3 non-spin-polarized x y
z '
endif
do i3=1,nr3
do i2=1,nr2
do i1=1,nr1
write(31,'(3i6,4(es24.14))')
i1,i2,i3,gridtt(i1,i2,i3),gridux(i1,i2,i3),griddy(i1,i2,i3),gridmz(i1,i2,i3)
end do
end do
end do
endif ! nspden
close(31)
exit
case(7)
open(unit=31,file=trim(filnam),form='unformatted')
xm=0 ; xp=rprimd(1,1)*Bohr_Ang
ym=0 ; yp=rprimd(2,2)*Bohr_Ang
zm=0 ; zp=rprimd(3,3)*Bohr_Ang
write(6, '(/,a,/)' )&
& ' Extremas (x,y,z) of the cube in which the molecule is placed, in
Angstroms'
write(6, '(5x,6f10.5)' ) xm,xp,ym,yp,zm,zp
write(6, '(/,a,2x,3i5)' )' Number of points per side: ',nr1,nr2,nr3
write(6, '(/,a,2x,i10,//)' )' Total number of points:', nr1*nr2*nr3
write(31) xm,xp,ym,yp,zm,zp,nr1,nr2,nr3
allocate(rhomacu(nr1,nr2))
do i3=1,nr3
do i2=1,nr2
do i1=1,nr1
rhomacu(i1,i2)=grid(i1,i2,i3)
end do
end do
write(31) rhomacu(:,:)
end do
close(31)
exit
case (8)
open(unit=31,file=trim(filnam),form='formatted')
write(6, '(/,a,/)' )&
& ' Extremas (x,y,z) of the cube in which the molecule is placed, in
Angstroms'
write(6, '(5x,6f10.5)' ) xm,xp,ym,yp,zm,zp
write(6, '(/,a,2x,3i5)' )' Number of points per side: ',nr1,nr2,nr3
write(6, '(/,a,2x,i10,//)' )' Total number of points:', nr1*nr2*nr3
write(31,*) 'TITLE = " " '
write(31,*) 'VARIABLES = "X" "Y" "Z" (all three in Angstrom)
"DENSITY or POTENTIAL" (atomic units) '
write(31,*) 'ZONE I=',nr1, ' J=', nr2, ' K=', nr3, ' F=POINT'
do i3=1,nr3
do i2=1,nr2
do i1=1,nr1
xnow = rprimd(1,1)*(i1-1)/nr1 + rprimd(1,2)*(i2-1)/nr2 +
rprimd(1,3)*(i3-1)/nr3
ynow = rprimd(2,1)*(i1-1)/nr1 + rprimd(2,2)*(i2-1)/nr2 +
rprimd(2,3)*(i3-1)/nr3
znow = rprimd(3,1)*(i1-1)/nr1 + rprimd(3,2)*(i2-1)/nr2 +
rprimd(3,3)*(i3-1)/nr3
write(31,*) Bohr_Ang*xnow, Bohr_Ang*ynow, Bohr_Ang*znow, grid
(i1,i2,i3)
end do
end do
end do
close(31)
exit
case (9)
open(unit=31,file=trim(filnam),form='formatted')
xm=0 ; xp=rprimd(1,1)*Bohr_Ang
ym=0 ; yp=rprimd(2,2)*Bohr_Ang
zm=0 ; zp=rprimd(3,3)*Bohr_Ang
write(6, '(/,a,/)' )&
& ' Extremas (x,y,z) of the cube in which the molecule is placed, in
Angstroms'
write(6, '(5x,6f10.5)' ) xm,xp,ym,yp,zm,zp
write(6, '(/,a,2x,3i5)' )' Number of points per side:
',nr1+1,nr2+1,nr3+1
write(6, '(/,a,2x,i10,//)' )' Total number of points:',
(nr1+1)*(nr2+1)*(nr3+1)
write(*,*) ' znucl = ', hdr%znucltypat, ' type = ', hdr%typat, '
ntypat = ', ntypat
gridshift1 = 0
gridshift2 = 0
gridshift3 = 0
write(*,*) 'Do you want to shift the grid along the x,y or z axis
(y/n)?'
write(*,*)
shift_tau(:) = zero
read (*,*) outputchar
if (outputchar == 'y' .or. outputchar == 'Y') then
write(*,*) 'Give the three shifts (x,y,z < ',nr1,nr2,nr3,') :'
write(*,*)
read (*,*) gridshift1, gridshift2, gridshift3
shift_tau(:) = gridshift1*rprimd(:,1)/(nr1+1) +
gridshift2*rprimd(:,2)/(nr2+1) + gridshift3*rprimd(:,3)/(nr3+1)
end if
!
! Generate translated coordinates to match density shift
!
allocate (tau2(3,natom))
do iatom = 1,natom
tau2(:,iatom) = xcart(:,iatom) - shift_tau(:)
end do
!DEBUG (LD)
!################################################################### (LD)
!TODO: document this bit: what is ispden 3???
!
if (ispden==3) then
! It is necessary to know previously how many atoms will be used.
! It is necessary to have an effective criteria in this point (not at
present),
! in order to plot the necessary magnetization arrows only.
index=natom
do i1=1,nr1
do i2=1,nr2
do i3=1,nr3
if(abs(gridux(i1,i2,i3)) .gt. 1.d-3 .and.&
abs(griddy(i1,i2,i3)) .gt. 1.d-3 .and. abs(gridmz(i1,i2,i3)) .gt.
1.d-3) then
index=index+1
endif
enddo
enddo
enddo
write(31,'(1X,A)') 'CRYSTAL'
write(31,'(1X,A)') 'PRIMVEC'
do i1 = 1,3
write(31,'(3(ES17.10,2X))') (Bohr_Ang*rprimd(i2,i1), i2=1,3)
end do
! write(31,'(1X,A)') 'BEGIN_BLOCK_DATAGRID3D'
! write(31,*) 'datagrids'
! write(31,'(1X,A)') 'DATAGRID_3D_DENSITY'
! write(31,*) nr1+1,nr2+1,nr3+1
! write(31,*) '0.0 0.0 0.0 '
! do i1 = 1,3
! write(31,'(3(ES17.10,2X))') (Bohr_Ang*rprimd(i2,i1), i2=1,3)
! end do
write(31,'(1X,A)') 'PRIMCOORD'
write(31,*) index, '1'
!
! write out atom types and positions
!
do iatom = 1,natom
write(31,'(i9,3(3X,ES17.10))') nint(hdr%znucltypat(hdr%typat(iatom))),&
& Bohr_Ang*tau2(1:3,iatom)
end do
! TODO: document
! are these magnetization vectors??
do i1=1,nr1
do i2=1,nr2
do i3=1,nr3
if(abs(gridux(i1,i2,i3)) .gt. 1.d-3 .and.&
abs(griddy(i1,i2,i3)) .gt. 1.d-3 .and. abs(gridmz(i1,i2,i3)) .gt.
1.d-3) then
write(31,'(A,1X,6(ES17.10,2X))')'X',i1/xp,i2/yp,i3/zp,&
gridux(i1,i2,i3),&
griddy(i1,i2,i3),&
gridmz(i1,i2,i3)
endif
enddo
enddo
enddo
! write (31,*)
! write(31,'(1X,A)') 'END_DATAGRID_3D'
! write(31,'(1X,A)') 'END_BLOCK_DATAGRID3D'
else
!################################################################### (LD)
!ENDDEBUG (LD)
!
! normal case: output density or potential (scalar field)
!
write(31,'(1X,A)') 'DIM-GROUP'
write(31,*) '3 1'
write(31,'(1X,A)') 'PRIMVEC'
do i1 = 1,3
write(31,'(3(ES17.10,2X))') (Bohr_Ang*rprimd(i2,i1), i2=1,3)
end do
write(31,'(1X,A)') 'PRIMCOORD'
write(31,*) natom, ' 1'
do iatom = 1,natom
write(31,'(i9,3(3X,ES17.10))') nint(hdr%znucltypat(hdr%typat(iatom))),&
& Bohr_Ang*tau2(1:3,iatom)
end do
write(31,'(1X,A)') 'ATOMS'
do iatom = 1,natom
write(31,'(i9,3(3X,ES17.10))') nint(hdr%znucltypat(hdr%typat(iatom))),&
& Bohr_Ang*tau2(1:3,iatom)
end do
! write(31,'(1X,A)') 'FRAMES'
write(31,'(1X,A)') 'BEGIN_BLOCK_DATAGRID3D'
write(31,*) 'datagrids'
write(31,'(1X,A)') 'DATAGRID_3D_DENSITY'
write(31,*) nr1+1,nr2+1,nr3+1
write(31,*) '0.0 0.0 0.0 '
do i1 = 1,3
write(31,'(3(ES17.10,2X))') (Bohr_Ang*rprimd(i2,i1), i2=1,3)
end do
index = 0
do ir3=gridshift3+1,nr3+1
ii3=mod(ir3-1,nr3) + 1
do ir2=gridshift2+1,nr2+1
ii2=mod(ir2-1,nr2) + 1
do ir1=gridshift1+1,nr1+1
ii1=mod(ir1-1,nr1) + 1
write(31,'(e20.10,2x)',ADVANCE='NO') grid(ii1,ii2,ii3)
index = index+1
if (mod (index,6) == 0) write (31,*)
end do
do ir1=1,gridshift1
ii1=mod(ir1-1,nr1) + 1
write(31,'(e20.10,2x)',ADVANCE='NO') grid(ii1,ii2,ii3)
index = index+1
if (mod (index,6) == 0) write (31,*)
end do
end do
do ir2=1,gridshift2
ii2=mod(ir2-1,nr2) + 1
do ir1=gridshift1+1,nr1+1
ii1=mod(ir1-1,nr1) + 1
write(31,'(e20.10,2x)',ADVANCE='NO') grid(ii1,ii2,ii3)
index = index+1
if (mod (index,6) == 0) write (31,*)
end do
do ir1=1,gridshift1
ii1=mod(ir1-1,nr1) + 1
write(31,'(e20.10,2x)',ADVANCE='NO') grid(ii1,ii2,ii3)
index = index+1
if (mod (index,6) == 0) write (31,*)
end do
end do
end do
do ir3=1,gridshift3
ii3=mod(ir3-1,nr3) + 1
do ir2=gridshift2+1,nr2+1
ii2=mod(ir2-1,nr2) + 1
do ir1=gridshift1+1,nr1+1
ii1=mod(ir1-1,nr1) + 1
write(31,'(e20.10,2x)',ADVANCE='NO') grid(ii1,ii2,ii3)
index = index+1
if (mod (index,6) == 0) write (31,*)
end do
do ir1=1,gridshift1
ii1=mod(ir1-1,nr1) + 1
write(31,'(e20.10,2x)',ADVANCE='NO') grid(ii1,ii2,ii3)
index = index+1
if (mod (index,6) == 0) write (31,*)
end do
end do
do ir2=1,gridshift2
ii2=mod(ir2-1,nr2) + 1
do ir1=gridshift1+1,nr1+1
ii1=mod(ir1-1,nr1) + 1
write(31,'(e20.10,2x)',ADVANCE='NO') grid(ii1,ii2,ii3)
index = index+1
if (mod (index,6) == 0) write (31,*)
end do
do ir1=1,gridshift1
ii1=mod(ir1-1,nr1) + 1
write(31,'(e20.10,2x)',ADVANCE='NO') grid(ii1,ii2,ii3)
index = index+1
if (mod (index,6) == 0) write (31,*)
end do
end do
end do
write (31,*)
write(31,'(1X,A)') 'END_DATAGRID_3D'
write(31,'(1X,A)') 'END_BLOCK_DATAGRID3D'
endif
close(31)
exit
case (10) ! formatted for OpenDX
open(unit=31,file=trim(filnam)//'.dx',form='formatted')
open(unit=32,file=trim(filnam)//'.xyz',form='formatted')
write(31, '(a,2x,3i5)' )'object 1 class gridpositions
counts',nr3,nr2,nr1
write(31, '(a)' )'origin 0 0 0'
write(31,'(a,3(es17.10,2x))')'delta ',(Bohr_Ang*rprimd(i1,3)/nr3,
i1=1,3)
write(31,'(a,3(es17.10,2x))')'delta ',(Bohr_Ang*rprimd(i1,2)/nr2,
i1=1,3)
write(31,'(a,3(es17.10,2x))')'delta ',(Bohr_Ang*rprimd(i1,1)/nr1,
i1=1,3)
write(31, '(a,2x,3i5)' )'object 2 class gridconnections
counts',nr3,nr2,nr1
write(31, '(a)' )'attribute "element type" string "cubes"'
write(31, '(a)' )'attribute "ref" string "positions"'
write(31, '(a,1x,i10,1x,a)' )'object 3 class array type float rank 0
items',nr1*nr2*nr3,' data follows'
do i3=1,nr3
do i2=1,nr2
do i1=1,nr1
write(31,'(es24.14)') grid(i1,i2,i3)
end do
end do
end do
write(31, '(a)' )'attribute "dep" string "positions"'
write(32, '(i6,/)' ) natom
do iatom=1,natom
do ii=1,3
xcart2(ii)=xcart(ii,iatom)
if (xred(ii,iatom)<-1e-4) then
xcart2(ii)=xcart2(ii)+rprimd(ii,ii)
else if (xred(ii,iatom)>rprimd(ii,ii)+1e-4) then
xcart2(ii)=xcart2(ii)-rprimd(ii,ii)
end if
end do
write(32,'(i8,3(3X,ES17.10))')
nint(hdr%znucltypat(hdr%typat(iatom))),Bohr_Ang*xcart2(1:3)
end do
write(31, '(a)' )'object "density" class field'
write(31, '(a)' )'component "positions" value 1'
write(31, '(a)' )'component "connections" value 2'
write(31, '(a)' )'component "data" value 3'
close(31)
close(32)
exit
case(11)
call
hirsh(grid,natom,nr1,nr2,nr3,ntypat,rprimd,xcart,hdr%typat,hdr%zionpsp,hdr%znucltypat)
exit
case(12)
#if defined HAVE_NETCDF
allocate(nuclz(natom))
filnam = trim(filnam)//'.nc'
! Creating the NetCDF file
status = nf90_create(filnam, nf90_clobber, ncid)
if (status /= nf90_noerr) call handle_ncerr(status, "Creating file")
! Ask for a title
write(*,*) 'Do you want a title in your NetCDF file? (0 = NO, 1 = YES)'
read(*,*) titlechoice
do
if (titlechoice ==0 .or. titlechoice ==1) exit
write(*,*) 'The answer is not correct, you must enter an integer
between 0 and 1 (0 = NO, 1 = YES)'
read(*,*) titlechoice
end do
if (titlechoice ==1) then
write(*,*) 'Enter your file''s title'
read(*,'(A)') filetitle
else
write(*,*) 'No title will be added in your NetCDF file'
end if
originatt(1:3,1:3)=0
do igrid = 1,3
gridwavefun1(igrid,1)=0
gridwavefun2(igrid,1)=0
gridwavefun3(igrid,1)=0
gridwavefun1(igrid,2)=Bohr_Ang*rprimd(igrid,3)/nr3
gridwavefun2(igrid,2)=Bohr_Ang*rprimd(igrid,2)/nr2
gridwavefun3(igrid,2)=Bohr_Ang*rprimd(igrid,1)/nr1
end do
do igrid =1,natom
nuclz(igrid)=hdr%znucltypat(hdr%typat(igrid))
end do
! Defining dimensions
status = nf90_def_dim(ncid,"gridsize1",nr1, gridsize1DimID)
if (status /= nf90_noerr) call handle_ncerr(status, "Defining
dimensions")
status = nf90_def_dim(ncid,"gridsize2",nr2, gridsize2DimID)
if (status /= nf90_noerr) call handle_ncerr(status, "Defining
dimensions")
status = nf90_def_dim(ncid,"gridsize3",nr3, gridsize3DimID)
if (status /= nf90_noerr) call handle_ncerr(status, "Defining
dimensions")
status = nf90_def_dim(ncid, "lat",3, latDimID)
if (status /= nf90_noerr) call handle_ncerr(status, "Defining
dimensions")
status = nf90_def_dim(ncid, "pos",2, posDimID)
if (status /= nf90_noerr) call handle_ncerr(status, "Defining
dimensions")
status = nf90_def_dim(ncid, "nbatom",natom, nbatomDimID)
if (status /= nf90_noerr) call handle_ncerr(status, "Defining
dimensions")
! Defining variables
status = nf90_def_var(ncid, "latticevec",nf90_float, (/ latDimID,
latDimID /), latticevecVarID)
if (status /= nf90_noerr) call handle_ncerr(status, "Defining
variables")
status = nf90_def_var(ncid, "origin",nf90_float, (/ latDimID, latDimID
/), originVarID)
if (status /= nf90_noerr) call handle_ncerr(status, "Defining
variables")
status = nf90_def_var(ncid, "atomposi",nf90_float, (/ latDimID,
nbatomDimID /), atomposiVarID)
if (status /= nf90_noerr) call handle_ncerr(status, "Defining
variables")
status = nf90_def_var(ncid, "atomicnum",nf90_float, nbatomDimID,
atomicnumVarID)
if (status /= nf90_noerr) call handle_ncerr(status, "Defining
variables")
status = nf90_def_var(ncid, "grid1",nf90_float, (/latDimID,posDimID/),
grid1VarID)
if (status /= nf90_noerr) call handle_ncerr(status, "Defining
variables")
status = nf90_def_var(ncid, "grid2",nf90_float, (/latDimID,posDimID/),
grid2VarID)
if (status /= nf90_noerr) call handle_ncerr(status, "Defining
variables")
status = nf90_def_var(ncid, "grid3",nf90_float, (/latDimID,posDimID/),
grid3VarID)
if (status /= nf90_noerr) call handle_ncerr(status, "Defining
variables")
if(fform0==51 .or. fform0==52)then ! Density case
status = nf90_def_var(ncid, "density",nf90_float, (/ gridsize1DimID,
gridsize2DimID, gridsize3DimID /), dataVarID)
if (status /= nf90_noerr) call handle_ncerr(status, "Defining
variables")
else if(fform0==101 .or. fform0==102)then ! Potential case
status =nf90_def_var(ncid, "potential",nf90_float, (/ gridsize1DimID,
gridsize2DimID, gridsize3DimID /), dataVarID)
if (status /= nf90_noerr) call handle_ncerr(status, "Defining
variables")
end if
! Defining attibutes
status = nf90_put_att(ncid,latticevecVarID , "field",
"latticevec,vector")
if (status /= nf90_noerr) call handle_ncerr(status, "Defining
attributes")
status = nf90_put_att(ncid,originVarID , "field", "origin,vector")
if (status /= nf90_noerr) call handle_ncerr(status, "Defining
attributes")
status = nf90_put_att(ncid,atomposiVarID , "field","atomposi,vector")
if (status /= nf90_noerr) call handle_ncerr(status, "Defining
attributes")
status = nf90_put_att(ncid,latticevecVarID , "positions","origin" )
if (status /= nf90_noerr) call handle_ncerr(status, "Defining
attributes")
status = nf90_put_att(ncid,atomicnumVarID , "positions", "atomposi")
if (status /= nf90_noerr) call handle_ncerr(status, "Defining
attributes")
status = nf90_put_att(ncid,dataVarID , "positions",
"grid1,product,compact;grid2,product,compact;grid3,product,compact")
if (status /= nf90_noerr) call handle_ncerr(status, "Defining
attributes")
status = nf90_put_att(ncid,dataVarID , "connections", (/nr3,nr2,nr1/))
if (status /= nf90_noerr) call handle_ncerr(status, "Defining
attributes")
status = nf90_put_att(ncid,latticevecVarID , "long_name", "Lattice
Vectors")
if (status /= nf90_noerr) call handle_ncerr(status, "Defining
attributes")
status = nf90_put_att(ncid,atomposiVarID , "long_name", "Atomic
Positions")
if (status /= nf90_noerr) call handle_ncerr(status, "Defining
attributes")
status = nf90_put_att(ncid, atomicnumVarID, "long_name", "Atomic
Numbers")
if (status /= nf90_noerr) call handle_ncerr(status, "Defining
attributes")
status = nf90_put_att(ncid,latticevecVarID, "units", "Angstroms")
if (status /= nf90_noerr) call handle_ncerr(status, "Defining
attributes")
status = nf90_put_att(ncid,atomposiVarID, "units", "Angstroms")
if (status /= nf90_noerr) call handle_ncerr(status, "Defining
attributes")
if(fform0==51 .or. fform0==52)then ! Density case
status = nf90_put_att(ncid, dataVarID, "long_name", "Density")
if (status /= nf90_noerr) call handle_ncerr(status, "Defining
attributes")
else if(fform0==101 .or. fform0==102)then ! Potential case
status = nf90_put_att(ncid, dataVarID, "long_name", "Potential")
if (status /= nf90_noerr) call handle_ncerr(status, "Defining
attributes")
end if
status = nf90_put_att(ncid, originVarID, "long_name", "Origin of the
Lattice")
if (status /= nf90_noerr) call handle_ncerr(status, "Defining
attributes")
if (titlechoice ==1) then
status = nf90_put_att(ncid, nf90_global, "title", filetitle)
if (status /= nf90_noerr) call handle_ncerr(status, "Defining
attributes")
end if
! Add the creation date
call date_and_time(strdat,strtime,strzone,values)
yyyy=values(1)
mm=values(2)
dd=values(3)
write(stridate(1:2),'(I2)') dd
stridate(3:3)=" "
stridate(4:7)=monnam(mm)
write(stridate(8:11),'(I4)') yyyy
status = nf90_put_att(ncid, nf90_global,"date", stridate)
if (status /= nf90_noerr) call handle_ncerr(status, "Defining
attributes")
! Ending the define mode and entering data mode
status = nf90_enddef(ncid)
if (status /= nf90_noerr) call handle_ncerr(status, "Entering data
mode")
! Putting the data
status = nf90_put_var(ncid,latticevecVarID,Bohr_Ang*rprimd)
if (status /= nf90_noerr) call handle_ncerr(status, "Putting data")
status = nf90_put_var(ncid,originVarID,originatt)
if (status /= nf90_noerr) call handle_ncerr(status, "Putting data")
status = nf90_put_var(ncid,atomposiVarID, Bohr_Ang*xcart)
if (status /= nf90_noerr) call handle_ncerr(status, "Putting data")
status = nf90_put_var(ncid,atomicnumVarID,nuclz)
if (status /= nf90_noerr) call handle_ncerr(status, "Putting data")
status = nf90_put_var(ncid,grid1VarID,gridwavefun1)
if (status /= nf90_noerr) call handle_ncerr(status, "Putting data")
status = nf90_put_var(ncid,grid2VarID,gridwavefun2)
if (status /= nf90_noerr) call handle_ncerr(status, "Putting data")
status = nf90_put_var(ncid,grid3VarID,gridwavefun3)
if (status /= nf90_noerr) call handle_ncerr(status, "Putting data")
status = nf90_put_var(ncid,dataVarID,grid)
if (status /= nf90_noerr) call handle_ncerr(status, "Putting data")
! Closing the file
status = nf90_close(ncid)
if (status /= nf90_noerr) call handle_ncerr(status, "Closing file")
write(*,*) 'The NetCDF file is done'
deallocate(nuclz)
write(*,*)
exit
#else
write(*,*) 'NetCDF is not defined. You must choose another option'
exit
#endif
case(0)
write(*,*)' Exit requested by user'
stop
case default
cycle
end select
end do
write(*,*) ' Task ',itask,' has been done !'
write(*,*)
write(*,*) ' More analysis of the 3D file ? (1=default=yes,0=no)'
read(*,*) iprompt
if(iprompt==0) then
exit
else
cycle
end if
end do
end if ! WF file or DEN/POT file
write(*,*)
write(*,*) ' Thank you for using me'
write(*,*)
!DEBUG
!write(*,*)' cut3d : before hdr_clean, densfileformat=',densfileformat
!write(*,*)' allocated(xcart)=',allocated(xcart)
!ENDDEBUG
if(densfileformat==1)then
call hdr_clean(hdr)
end if
deallocate(xcart)
end program cut3d
!!***
- Re: [abinit-forum] Problem with cut3d and Xcrysden format, Matthieu Verstraete, 07/06/2007
- Re: [abinit-forum] Problem with cut3d and Xcrysden format, Mauro Sgroi, 07/06/2007
Archive powered by MHonArc 2.6.16.