scnc
/
quantum-espresso
mirror of https://gitlab.com/QEF/q-e.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity
quantum-espresso/CPV/fpmdpp.f90

962 lines
29 KiB
Fortran
Raw Blame History

!
! Copyright (C) 2002-2008 Quantum ESPRESSO group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!
! This file holds XSF (=Xcrysden Structure File) utilities.
! Routines written by Tone Kokalj on Mon Jan 27 18:51:17 CET 2003
! modified by Gerardo Ballabio and Carlo Cavazzoni
! on Thu Jul 22 18:57:26 CEST 2004
!
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
! --------------------------------------------------------------------
! this routine writes the crystal structure in XSF, GRD and PDB format
! from a FPMD output files
! --------------------------------------------------------------------
PROGRAM fpmd_postproc
USE kinds, ONLY : DP
USE constants, ONLY : bohr => BOHR_RADIUS_ANGS
USE io_files, ONLY : tmp_dir, prefix, iunpun, xmlpun, outdir
USE io_global, ONLY : io_global_start, io_global_getmeta
USE mp_global, ONLY : mp_global_start, init_pool
USE mp, ONLY : mp_end, mp_start, mp_env
USE iotk_module
USE xml_io_base
IMPLICIT NONE
INTEGER, PARAMETER :: maxsp = 20
INTEGER :: natoms, nsp, na(maxsp), atomic_number(maxsp)
INTEGER :: ounit, cunit, punit, funit, dunit, bunit, ksunit
INTEGER :: nr1, nr2, nr3, ns1, ns2, ns3
INTEGER :: np1, np2, np3, np, ispin
INTEGER, ALLOCATABLE :: ityp(:)
REAL(DP) :: at(3, 3), atinv(3, 3), ht0(3, 3), h0(3, 3)
REAL(DP) :: rhof, rhomax, rhomin, rhoc(6)
REAL(DP), ALLOCATABLE :: rho_in(:,:,:), rho_out(:,:,:)
REAL(DP), ALLOCATABLE :: tau_in(:,:), tau_out(:,:)
REAL(DP), ALLOCATABLE :: sigma(:,:), force(:,:)
REAL(DP), ALLOCATABLE :: stau0(:,:), svel0(:,:), force0(:,:)
CHARACTER(len=256) :: filepp, fileout, output
CHARACTER(len=256) :: filecel, filepos, filefor, filepdb
CHARACTER(len=256) :: print_state
CHARACTER(len=3) :: atm( maxsp ), lab
CHARACTER(len=4) :: charge_density
LOGICAL :: lcharge, lforces, ldynamics, lpdb, lrotation
LOGICAL :: lbinary, found
INTEGER :: nframes
INTEGER :: ios, nat, ndr
INTEGER :: nproc, mpime, world, root
REAL(DP) :: x, y, z, fx, fy, fz
INTEGER :: i, j, k, n, ix, iy, iz, ierr
REAL(DP) :: euler(6)
NAMELIST /inputpp/ prefix, fileout, output, outdir, &
lcharge, lforces, ldynamics, lpdb, lrotation, &
ns1, ns2, ns3, np1, np2, np3, print_state, &
atomic_number, nframes, ndr, charge_density, &
lbinary
! default values
dunit = 14
! ... Intel compilers v .ge.8 allocate a lot of stack space
! ... Stack limit is often small, thus causing SIGSEGV and crash
CALL remove_stack_limit ( )
! see cprstart.f90 for the meaning of the following 4 calls
CALL mp_start()
CALL mp_env( nproc, mpime, world )
CALL mp_global_start( root, mpime, world, nproc )
CALL io_global_start( mpime, root )
CALL get_env( 'ESPRESSO_TMPDIR', outdir )
IF ( TRIM( outdir ) == ' ' ) outdir = './'
prefix = 'cp'
fileout = 'out'
output = 'xsf' ! 'grd'
outdir = './'
lcharge = .false.
lforces = .false.
ldynamics = .false.
lpdb = .false.
lrotation = .false.
ns1 = 0
ns2 = 0
ns3 = 0
np1 = 1
np2 = 1
np3 = 1 !
nframes = 1 ! number of MD step to be read to buind the trajectory
ndr = 51 ! restart file number
atomic_number = 1 ! atomic number of the species in the restart file
charge_density = 'full' ! specify the component to plot: 'full', 'spin'
print_state = ' ' ! specify the Kohn-Sham state to plot: 'KS_1'
lbinary = .TRUE.
call input_from_file()
! read namelist
READ( 5, inputpp, iostat=ios)
! set file names
!
filecel = TRIM(outdir) // TRIM(prefix) // '.cel'
filepos = TRIM(outdir) // TRIM(prefix) // '.pos'
filefor = TRIM(outdir) // TRIM(prefix) // '.for'
!
filepdb = TRIM(fileout) // '.pdb'
!
! append extension
!
IF (output == 'xsf') THEN
IF (ldynamics) THEN
fileout = TRIM(fileout) // '.axsf'
ELSE
fileout = TRIM(fileout) // '.xsf'
END IF
ELSE IF (output == 'grd') THEN
fileout = TRIM(fileout) // '.grd'
END IF
np = np1 * np2 * np3
IF (np1 < 1 .OR. np2 < 1 .OR. np3 < 1) THEN
WRITE(*,*) 'Error: zero or negative replicas not allowed'
STOP
END IF
! check for wrong input
IF (ldynamics .AND. nframes < 2) THEN
WRITE(*,*) 'Error: dynamics requested, but only one frame'
STOP
END IF
IF (.NOT. ldynamics) nframes = 1
IF (ldynamics .AND. lcharge) THEN
WRITE(*,*) 'Error: dynamics with charge density not supported'
STOP
END IF
IF (ldynamics .AND. ( print_state /= ' ' ) ) THEN
WRITE(*,*) 'Error: dynamics with print_state not supported'
STOP
END IF
!
! Now read the XML data file
!
filepp = restart_dir( outdir, ndr )
!
filepp = TRIM( filepp ) // '/' // TRIM(xmlpun)
!
CALL iotk_open_read( dunit, file = TRIM( filepp ), BINARY = .FALSE., ROOT = attr )
CALL iotk_scan_begin( dunit, "IONS", FOUND = found )
IF( .NOT. found ) THEN
CALL errore( ' cppp ', ' IONS not found in data-file.xml ', 1 )
END IF
CALL iotk_scan_dat( dunit, "NUMBER_OF_ATOMS", nat )
CALL iotk_scan_dat( dunit, "NUMBER_OF_SPECIES", nsp )
ALLOCATE( ityp( nat * np ) ) ! atomic species
DO i = 1, nsp
!
CALL iotk_scan_begin( dunit, "SPECIE" // TRIM( iotk_index( i ) ), FOUND = found )
!
IF( .NOT. found ) THEN
CALL errore( ' cppp ', "SPECIE" // TRIM( iotk_index( i ) ) // ' not found in data-file.xml ', 1 )
END IF
CALL iotk_scan_dat( dunit, "ATOM_TYPE", atm(i) )
!
! CALL iotk_scan_dat( dunit, &
! TRIM( atm(i) )//"_MASS", amass(i), ATTR = attr )
!
! CALL iotk_scan_dat( dunit, &
! "PSEUDO_FOR_" // TRIM( atm(i) ), psfile(i) )
!
CALL iotk_scan_end( dunit, "SPECIE" // TRIM( iotk_index( i ) ) )
!
END DO
!
! CALL iotk_scan_empty( dunit, "UNITS_FOR_ATOMIC_POSITIONS", attr )
! CALL iotk_scan_attr( attr, "UNIT", pos_unit )
!
DO i = 1, nat
!
CALL iotk_scan_empty( dunit, "ATOM" // TRIM( iotk_index( i ) ), attr )
CALL iotk_scan_attr( attr, "SPECIES", lab )
CALL iotk_scan_attr( attr, "INDEX", ityp(i) )
! CALL iotk_scan_attr( attr, "tau", tau(:,i) )
! CALL iotk_scan_attr( attr, "if_pos", if_pos(:,i) )
!
END DO
CALL iotk_scan_end( dunit, "IONS" )
CALL iotk_scan_begin( dunit, "PLANE_WAVES" )
CALL iotk_scan_empty( dunit, "FFT_GRID", attr )
CALL iotk_scan_attr( attr, "nr1", nr1 )
CALL iotk_scan_attr( attr, "nr2", nr2 )
CALL iotk_scan_attr( attr, "nr3", nr3 )
CALL iotk_scan_end( dunit, "PLANE_WAVES" )
ALLOCATE( stau0( 3, nat ) )
ALLOCATE( svel0( 3, nat ) )
ALLOCATE( force0( 3, nat ) ) ! forces, atomic units
CALL iotk_scan_begin( dunit, "TIMESTEPS", attr )
CALL iotk_scan_begin( dunit, "STEP0" )
CALL iotk_scan_begin( dunit, "IONS_POSITIONS" )
CALL iotk_scan_dat( dunit, "stau", stau0 )
CALL iotk_scan_dat( dunit, "svel", svel0 )
CALL iotk_scan_dat( dunit, "force", force0 )
CALL iotk_scan_end( dunit, "IONS_POSITIONS" )
CALL iotk_scan_begin( dunit, "CELL_PARAMETERS" )
CALL iotk_scan_dat( dunit, "ht", ht0 )
CALL iotk_scan_end( dunit, "CELL_PARAMETERS" )
CALL iotk_scan_end( dunit, "STEP0" )
CALL iotk_scan_end( dunit, "TIMESTEPS" )
!
ispin = 1
!
!
CALL iotk_close_read( dunit )
!
! End of reading from data file
!
IF ( nsp > maxsp ) THEN
WRITE(*,*) 'Error: too many atomic species'
STOP
END IF
natoms = nat
!
! Count atoms in each species
!
na = 0
DO i = 1, nat
na( ityp( i ) ) = na( ityp( i ) ) + 1 ! total number of atoms
END DO
! assign species (from input) to each atom
!
k = 0
DO i = 1, nsp
DO j = 1, na(i)
k = k + 1
ityp(k) = atomic_number(i)
END DO
END DO
! allocate arrays
ALLOCATE(tau_in(3, nat)) ! atomic positions, angstroms
ALLOCATE(tau_out(3, nat * np)) ! replicated positions
ALLOCATE(sigma(3, nat ) ) ! scaled coordinates
!
IF (lforces) ALLOCATE( force( 3, nat * np ) )
! charge density
IF ( lcharge .OR. print_state /= ' ' ) THEN
IF (ns1 == 0) ns1 = nr1
IF (ns2 == 0) ns2 = nr2
IF (ns3 == 0) ns3 = nr3
ALLOCATE( rho_in ( nr1, nr2, nr3 ) ) ! original charge density
ALLOCATE( rho_out( ns1, ns2, ns3 ) ) ! rescaled charge density
END IF
! open output file for trajectories or charge density
!
ounit = 10
OPEN(ounit, file=fileout, status='unknown')
! open Cell trajectory file
!
cunit = 11
OPEN(cunit, file=filecel, status='old')
! open Positions trajectory file
!
punit = 12
OPEN(punit, file=filepos, status='old')
! open Force trajectory file
!
funit = 13
if (lforces) OPEN(funit, file=filefor, status='old')
! open PDB file
!
bunit = 15
OPEN(bunit, file=filepdb, status='unknown')
! Unit for KS states
!
ksunit = 16
! XSF file header
!
IF ( output == 'xsf' ) THEN
IF ( ldynamics ) WRITE(ounit,*) 'ANIMSTEPS', nframes
WRITE( ounit, * ) 'CRYSTAL'
END IF
DO n = 1, nframes
!
IF ( ldynamics ) WRITE(*,'("frame",1X,I4)') n
! read data from files produced by fpmd
!
CALL read_fpmd( lforces, lcharge, lbinary, cunit, punit, funit, dunit, &
natoms, nr1, nr2, nr3, ispin, at, tau_in, force, &
rho_in, prefix, outdir, ndr, charge_density )
IF( nframes == 1 ) THEN
!
! use values from the XML file
!
IF( lforces ) force( 1:3, 1:nat ) = force0( 1:3, 1:nat )
!
h0 = TRANSPOSE( ht0 )
!
! from scaled to real coordinates
!
tau_in( :, : ) = MATMUL( h0( :, : ), stau0( :, : ) )
!
! convert atomic units to Angstroms
!
at = h0 * bohr
tau_in = tau_in * bohr
!
END IF
WRITE(*,'(2x,"Cell parameters (Angstroms):")')
WRITE(*,'(3(2x,f10.6))') ((at(i, j), i=1,3), j=1,3)
!
WRITE(*,'(2x,"Atomic coordinates (Angstroms):")')
WRITE(*,'(3(2x,f10.6))') ((tau_in(i, j), i=1,3), j=1,natoms)
! compute scaled coordinates
!
CALL inverse( at, atinv )
sigma(:,:) = MATMUL(atinv(:,:), tau_in(:,:))
! compute cell dimensions and Euler angles
CALL at_to_euler( at, euler )
IF (lpdb) THEN
! apply periodic boundary conditions
DO i = 1, natoms
DO j = 1, 3
sigma(j, i) = sigma(j, i) - FLOOR(sigma(j, i))
END DO
END DO
! recompute Cartesian coordinates
tau_in(:,:) = MATMUL(at(:,:), sigma(:,:))
END IF
IF (lrotation) THEN
! compute rotated cell
CALL euler_to_at( euler, at )
! rotate atomic positions as well
tau_in(:,:) = MATMUL(at(:,:), sigma(:,:))
END IF
! replicate atoms
k = 0
DO ix = 1, np1
DO iy = 1, np2
DO iz = 1, np3
DO j = 1, natoms
k = k + 1
tau_out(:, k) = tau_in(:, j) + (ix-1) * at(:, 1) + &
(iy-1) * at(:, 2) + (iz-1) * at(:, 3)
ityp(k) = ityp(j)
IF (lforces) force(:, k) = force(:, j)
END DO
END DO
END DO
END DO
natoms = natoms * np
! compute supercell
at(:, 1) = at(:, 1) * np1
at(:, 2) = at(:, 2) * np2
at(:, 3) = at(:, 3) * np3
euler(1) = euler(1) * np1
euler(2) = euler(2) * np2
euler(3) = euler(3) * np3
IF ( lcharge ) &
CALL scale_charge( rho_in, rho_out, nr1, nr2, nr3, ns1, ns2, ns3, &
np1, np2, np3 )
IF ( output == 'xsf' ) THEN
! write data as XSF format
CALL write_xsf( ldynamics, lforces, lcharge, ounit, n, at, &
natoms, ityp, tau_out, force, rho_out, &
ns1, ns2, ns3 )
ELSE IF( output == 'grd' ) THEN
! write data as GRD format
CALL write_grd( ounit, at, rho_out, ns1, ns2, ns3 )
END IF
END DO
CLOSE(ounit)
IF ( print_state /= ' ' ) THEN
!
CALL read_density( TRIM( print_state ) // '.xml', dunit, nr1, nr2, nr3, rho_in, lbinary )
CALL scale_charge( rho_in, rho_out, nr1, nr2, nr3, ns1, ns2, ns3, np1, np2, np3 )
!
IF (output == 'xsf') THEN
! write data as XSF format
OPEN( unit = ksunit, file = TRIM( print_state ) // '.xsf' )
WRITE( ksunit, * ) 'CRYSTAL' ! XSF files need this one line header
CALL write_xsf( ldynamics, lforces, .true., ksunit, n, at, &
natoms, ityp, tau_out, force, rho_out, ns1, ns2, ns3 )
ELSE IF( output == 'grd' ) THEN
OPEN( unit = ksunit, file = TRIM( print_state ) // '.grd' )
CALL write_grd( ksunit, at, rho_out, ns1, ns2, ns3 )
END IF
!
CLOSE( ksunit )
!
END IF
! write atomic positions as PDB format
CALL write_pdb( bunit, at, tau_out, natoms, ityp, euler, lrotation )
! free allocated resources
CLOSE(punit)
CLOSE(cunit)
IF (lforces) CLOSE(funit)
DEALLOCATE(tau_in)
DEALLOCATE(tau_out)
DEALLOCATE(ityp)
IF( ALLOCATED( force ) ) DEALLOCATE(force)
IF( ALLOCATED( rho_in ) ) DEALLOCATE(rho_in)
IF( ALLOCATED( rho_out ) ) DEALLOCATE(rho_out)
DEALLOCATE( stau0 )
DEALLOCATE( svel0 )
DEALLOCATE( force0 )
CALL mp_end()
STOP
END PROGRAM fpmd_postproc
!
!
!
SUBROUTINE read_fpmd( lforces, lcharge, lbinary, cunit, punit, funit, dunit, &
natoms, nr1, nr2, nr3, ispin, at, tau, force, &
rho, prefix, outdir, ndr, charge_density )
USE kinds, ONLY: DP
USE constants, ONLY: bohr => BOHR_RADIUS_ANGS
USE xml_io_base
USE iotk_module
IMPLICIT NONE
LOGICAL, INTENT(in) :: lforces, lcharge, lbinary
INTEGER, INTENT(in) :: cunit, punit, funit, dunit
INTEGER, INTENT(in) :: natoms, nr1, nr2, nr3, ispin, ndr
REAL(DP), INTENT(out) :: at(3, 3), tau(3, natoms), force(3, natoms)
REAL(DP), INTENT(out) :: rho(nr1, nr2, nr3)
CHARACTER(LEN=*), INTENT(IN) :: prefix
CHARACTER(LEN=*), INTENT(IN) :: outdir
CHARACTER(LEN=*), INTENT(IN) :: charge_density
INTEGER :: i, j, ix, iy, iz, ierr
REAL(DP) :: rhomin, rhomax, rhof
REAL(DP) :: x, y, z, fx, fy, fz
CHARACTER(LEN=256) :: filename
INTEGER :: n1, n2, n3
REAL(DP), ALLOCATABLE :: rho_plane(:)
! read cell vectors
! NOTE: colums are lattice vectors
!
READ(cunit,*)
DO i = 1, 3
READ(cunit,*) ( at(i, j), j=1,3 )
END DO
at(:, :) = at(:, :) * bohr
! read atomic coordinates
READ(punit,*)
IF (lforces) READ(funit,*)
DO i = 1, natoms
! convert atomic units to Angstroms
READ(punit,*) x, y, z
tau(1, i) = x * bohr
tau(2, i) = y * bohr
tau(3, i) = z * bohr
IF (lforces) THEN
! read forces
READ (funit,*) fx, fy, fz
force(1, i) = fx
force(2, i) = fy
force(3, i) = fz
END IF
END DO
IF (lcharge) THEN
filename = restart_dir( outdir, ndr )
!
IF( charge_density == 'spin' ) THEN
filename = TRIM( filename ) // '/' // 'spin-polarization'
ELSE
filename = TRIM( filename ) // '/' // 'charge-density'
END IF
!
!
IF ( check_file_exst ( TRIM(filename)//'.dat' ) ) THEN
!
CALL read_density( TRIM(filename)//'.dat', dunit, nr1, nr2, nr3, rho, lbinary )
!
ELSEIF ( check_file_exst ( TRIM(filename)//'.xml' ) ) THEN
!
CALL read_density( TRIM(filename)//'.xml', dunit, nr1, nr2, nr3, rho, lbinary )
!
ELSE
CALL infomsg ('read_fpmd', 'file '//TRIM(filename)//' not found' )
ENDIF
!
END IF
RETURN
END SUBROUTINE read_fpmd
SUBROUTINE read_density( filename, dunit, nr1, nr2, nr3, rho, lbinary )
USE kinds, ONLY: DP
USE xml_io_base
USE iotk_module
IMPLICIT NONE
LOGICAL, INTENT(in) :: lbinary
INTEGER, INTENT(in) :: dunit
INTEGER, INTENT(in) :: nr1, nr2, nr3
REAL(DP), INTENT(out) :: rho(nr1, nr2, nr3)
CHARACTER(LEN=*), INTENT(IN) :: filename
INTEGER :: ix, iy, iz, ierr
REAL(DP) :: rhomin, rhomax, rhof
INTEGER :: n1, n2, n3
REAL(DP), ALLOCATABLE :: rho_plane(:)
!
WRITE(*,'("Reading density from: ", A80)' ) TRIM( filename )
!
CALL iotk_open_read( dunit, file = TRIM( filename ) , BINARY = lbinary, ROOT = attr, IERR = ierr )
!
CALL iotk_scan_begin( dunit, "CHARGE-DENSITY" )
CALL iotk_scan_empty( dunit, "INFO", attr )
CALL iotk_scan_attr( attr, "nr1", n1 )
CALL iotk_scan_attr( attr, "nr2", n2 )
CALL iotk_scan_attr( attr, "nr3", n3 )
!
ALLOCATE( rho_plane( n1 * n2 ) )
! read charge density from file
! note: must transpose
DO iz = 1, n3
CALL iotk_scan_dat( dunit, "z" // iotk_index( iz ), rho_plane )
IF( iz <= nr3 ) THEN
DO iy = 1, MIN( n2, nr2 )
DO ix = 1, MIN( n1, nr1 )
rho(ix, iy, iz) = rho_plane( ix + ( iy - 1 ) * n1 )
END DO
END DO
END IF
END DO
CALL iotk_scan_end( dunit, "CHARGE-DENSITY" )
CALL iotk_close_read( dunit )
rhomin = MINVAL(rho(:,:,:))
rhomax = MAXVAL(rho(:,:,:))
! print some info
WRITE(*,'(2x,"Density grid:")')
WRITE(*,'(3(2x,i6))') nr1, nr2, nr3
WRITE(*,'(2x,"spin = ",A4)') filename
WRITE(*,'(2x,"Minimum and maximum values:")')
WRITE(*,'(3(2x,1pe12.4))') rhomin, rhomax
RETURN
END SUBROUTINE read_density
!
!
!
! compute inverse of 3*3 matrix
!
SUBROUTINE inverse( at, atinv )
IMPLICIT NONE
INTEGER, PARAMETER :: DP = KIND(0.0d0)
REAL(DP), INTENT(in) :: at(3, 3)
REAL(DP), INTENT(out) :: atinv(3, 3)
REAL(DP) :: det
atinv(1, 1) = at(2, 2) * at(3, 3) - at(2, 3) * at(3, 2)
atinv(2, 1) = at(2, 3) * at(3, 1) - at(2, 1) * at(3, 3)
atinv(3, 1) = at(2, 1) * at(3, 2) - at(2, 2) * at(3, 1)
atinv(1, 2) = at(1, 3) * at(3, 2) - at(1, 2) * at(3, 3)
atinv(2, 2) = at(1, 1) * at(3, 3) - at(1, 3) * at(3, 1)
atinv(3, 2) = at(1, 2) * at(3, 1) - at(1, 1) * at(3, 2)
atinv(1, 3) = at(1, 2) * at(2, 3) - at(1, 3) * at(2, 2)
atinv(2, 3) = at(1, 3) * at(2, 1) - at(1, 1) * at(2, 3)
atinv(3, 3) = at(1, 1) * at(2, 2) - at(1, 2) * at(2, 1)
det = at(1, 1) * atinv(1, 1) + at(1, 2) * atinv(2, 1) + &
at(1, 3) * atinv(3, 1)
atinv(:,:) = atinv(:,:) / det;
RETURN
END SUBROUTINE inverse
! generate cell dimensions and Euler angles from cell vectors
! euler(1:6) = a, b, c, alpha, beta, gamma
! I didn't call the array "celldm" because that could be confusing,
! since in PWscf the convention is different:
! celldm(1:6) = a, b/a, c/a, cos(alpha), cos(beta), cos(gamma)
SUBROUTINE at_to_euler( at, euler )
IMPLICIT NONE
INTEGER, PARAMETER :: DP = KIND(0.0d0)
REAL(DP), INTENT(in) :: at(3, 3)
REAL(DP), INTENT(out) :: euler(6)
REAL(DP), PARAMETER :: rad2deg = 180.0d0 / 3.14159265358979323846d0
REAL(DP) :: dot(3, 3)
INTEGER :: i, j
DO i = 1, 3
DO j = i, 3
dot(i, j) = dot_product(at(:,i), at(:,j))
END DO
END DO
DO i = 1, 3
euler(i) = sqrt(dot(i, i))
END DO
euler(4) = acos(dot(2, 3) / (euler(2) * euler(3))) * rad2deg
euler(5) = acos(dot(1, 3) / (euler(1) * euler(3))) * rad2deg
euler(6) = acos(dot(1, 2) / (euler(1) * euler(2))) * rad2deg
RETURN
END SUBROUTINE at_to_euler
! generate cell vectors back from cell dimensions and Euler angles
! euler(1:6) = a, b, c, alpha, beta, gamma
! here I follow the PDB convention, namely, c is oriented along the z
! axis and b lies in the yz plane, or to put it another way, at is
! lower triangular
SUBROUTINE euler_to_at( euler, at )
IMPLICIT NONE
INTEGER, PARAMETER :: DP = KIND(0.0d0)
REAL(DP), PARAMETER :: deg2rad = 3.14159265358979323846d0 / 180.0d0
REAL(DP), INTENT(in) :: euler(6)
REAL(DP), INTENT(out) :: at(3, 3)
REAL(DP) :: cos_ab, cos_ac, cos_bc, temp1, temp2
cos_bc = COS(euler(4) * deg2rad)
cos_ac = COS(euler(5) * deg2rad)
cos_ab = COS(euler(6) * deg2rad)
temp1 = SQRT(1.0d0 - cos_bc*cos_bc) ! sin_bc
temp2 = (cos_ab - cos_bc*cos_ac) / temp1
at(1, 1) = SQRT(1.0d0 - cos_ac*cos_ac - temp2*temp2) * euler(1)
at(2, 1) = temp2 * euler(1)
at(3, 1) = cos_ac * euler(1)
at(1, 3) = 0.0d0
at(2, 3) = 0.0d0
at(3, 3) = euler(3)
at(1, 2) = 0.0d0
at(2, 2) = temp1 * euler(2)
at(3, 2) = cos_bc * euler(2)
RETURN
END SUBROUTINE euler_to_at
! map charge density from a grid to another by linear interpolation
! along the three axes
SUBROUTINE scale_charge( rho_in, rho_out, nr1, nr2, nr3, ns1, ns2, ns3, &
np1, np2, np3 )
IMPLICIT NONE
INTEGER, PARAMETER :: DP = KIND(0.0d0)
INTEGER, INTENT(in) :: nr1, nr2, nr3, ns1, ns2, ns3, np1, np2, np3
REAL(DP), INTENT(in) :: rho_in( nr1, nr2, nr3 )
REAL(DP), INTENT(out) :: rho_out( ns1, ns2, ns3 )
INTEGER :: i, j, k
INTEGER :: i0(ns1), j0(ns2), k0(ns3), i1(ns1), j1(ns2), k1(ns3)
REAL(DP) :: x0(ns1), y0(ns2), z0(ns3), x1(ns1), y1(ns2), z1(ns3)
! precompute interpolation data
DO i = 1, ns1
CALL scale_linear( i, nr1, ns1, np1, i0(i), i1(i), x0(i), x1(i) )
END DO
DO j = 1, ns2
CALL scale_linear( j, nr2, ns2, np2, j0(j), j1(j), y0(j), y1(j) )
END DO
DO k = 1, ns3
CALL scale_linear( k, nr3, ns3, np3, k0(k), k1(k), z0(k), z1(k) )
END DO
! interpolate linearly along three axes
DO i = 1, ns1
DO j = 1, ns2
DO k = 1, ns3
rho_out(i, j, k) = &
rho_in(i1(i), j1(j), k1(k)) * x0(i) * y0(j) * z0(k) + &
rho_in(i0(i), j1(j), k1(k)) * x1(i) * y0(j) * z0(k) + &
rho_in(i1(i), j0(j), k1(k)) * x0(i) * y1(j) * z0(k) + &
rho_in(i1(i), j1(j), k0(k)) * x0(i) * y0(j) * z1(k) + &
rho_in(i0(i), j0(j), k1(k)) * x1(i) * y1(j) * z0(k) + &
rho_in(i0(i), j1(j), k0(k)) * x1(i) * y0(j) * z1(k) + &
rho_in(i1(i), j0(j), k0(k)) * x0(i) * y1(j) * z1(k) + &
rho_in(i0(i), j0(j), k0(k)) * x1(i) * y1(j) * z1(k)
END DO
END DO
END DO
RETURN
END SUBROUTINE scale_charge
! compute grid parameters for linear interpolation
SUBROUTINE scale_linear( n, nr, ns, np, n0, n1, r0, r1 )
IMPLICIT NONE
INTEGER, PARAMETER :: DP = KIND(0.0d0)
INTEGER, INTENT(in) :: n, nr, ns, np
INTEGER, INTENT(out) :: n0, n1
REAL(DP), INTENT(out) :: r0, r1
! map new grid point onto old grid
! mapping is: 1 --> 1, ns+1 --> (nr*np)+1
r0 = REAL((n-1) * nr*np, DP) / ns + 1.0d0
! indices of neighbors
n0 = int(r0)
n1 = n0 + 1
! distances from neighbors
r0 = r0 - n0
r1 = 1.0d0 - r0
! apply periodic boundary conditions
n0 = MOD(n0 - 1, nr) + 1
n1 = MOD(n1 - 1, nr) + 1
RETURN
END SUBROUTINE scale_linear
SUBROUTINE write_xsf( ldynamics, lforces, lcharge, ounit, n, at, &
natoms, ityp, tau, force, rho, nr1, nr2, nr3 )
IMPLICIT NONE
INTEGER, PARAMETER :: DP = KIND(0.0d0)
LOGICAL, INTENT(in) :: ldynamics, lforces, lcharge
INTEGER, INTENT(in) :: ounit, n, natoms, ityp(natoms)
INTEGER, INTENT(in) :: nr1, nr2, nr3
REAL(DP), INTENT(in) :: at(3, 3), tau(3, natoms), force(3, natoms)
REAL(DP), INTENT(in) :: rho(nr1, nr2, nr3)
INTEGER :: i, j, ix, iy, iz
! write cell
IF (ldynamics) THEN
WRITE(ounit,*) 'PRIMVEC', n
ELSE
WRITE(ounit,*) 'PRIMVEC'
END IF
WRITE(ounit,'(2(3f15.9/),3f15.9)') at
IF (ldynamics) THEN
WRITE(ounit,*) 'CONVVEC', n
WRITE(ounit,'(2(3f15.9/),3f15.9)') at
END IF
! write atomic coordinates (and forces)
IF (ldynamics) THEN
WRITE(ounit,*) 'PRIMCOORD', n
ELSE
WRITE(ounit,*) 'PRIMCOORD'
END IF
WRITE(ounit,*) natoms, 1
DO i = 1, natoms
IF (lforces) THEN
WRITE (ounit,'(i3,3x,3f15.9,1x,3f12.5)') ityp(i), &
(tau(j, i), j=1,3), (force(j, i), j=1,3)
ELSE
WRITE (ounit,'(i3,3x,3f15.9,1x,3f12.5)') ityp(i), &
(tau(j, i), j=1,3)
END IF
END DO
! write charge density
IF (lcharge) THEN
! XSF scalar-field header
WRITE(ounit,'(a)') 'BEGIN_BLOCK_DATAGRID_3D'
WRITE(ounit,'(a)') '3D_PWSCF'
WRITE(ounit,'(a)') 'DATAGRID_3D_UNKNOWN'
! mesh dimensions
WRITE(ounit,*) nr1, nr2, nr3
! origin
WRITE(ounit,'(3f10.6)') 0.0d0, 0.0d0, 0.0d0
! lattice vectors
WRITE(ounit,'(3f10.6)') ((at(i, j), i=1,3), j=1,3)
! charge density
WRITE(ounit,'(6e13.5)') &
(((rho(ix, iy, iz), ix=1,nr1), iy=1,nr2), iz=1,nr3)
WRITE(ounit,'(a)') 'END_DATAGRID_3D'
WRITE(ounit,'(a)') 'END_BLOCK_DATAGRID_3D'
END IF
RETURN
END SUBROUTINE write_xsf
SUBROUTINE write_grd( ounit, at, rho, nr1, nr2, nr3 )
IMPLICIT NONE
INTEGER, PARAMETER :: DP = KIND(0.0d0)
INTEGER, INTENT(in) :: ounit
INTEGER, INTENT(in) :: nr1, nr2, nr3
REAL(DP), INTENT(in) :: at(3, 3), rho(nr1, nr2, nr3)
INTEGER :: i, j, k
REAL(DP) :: euler(6)
CALL at_to_euler( at, euler )
WRITE(ounit,*) 'charge density'
WRITE(ounit,*) '(1p,e12.5)'
WRITE(ounit,fmt='(6f9.3)') (euler(i), i=1,6)
WRITE(ounit,fmt='(3i5)') nr1 - 1, nr2 - 1, nr3 - 1
WRITE(ounit,fmt='(7i5)') 1, 0, 0, 0, nr1 - 1, nr2 - 1, nr3 - 1
WRITE(ounit,fmt='(1p,e12.5)') (((rho(i, j, k), i=1,nr1), j=1,nr2), k=1,nr3)
RETURN
END SUBROUTINE write_grd
SUBROUTINE write_pdb( bunit, at, tau, natoms, ityp, euler, lrotation )
IMPLICIT NONE
INTEGER, PARAMETER :: DP = KIND(0.0d0)
INTEGER, INTENT(in) :: bunit, natoms
INTEGER, INTENT(in) :: ityp(natoms)
REAL(DP), INTENT(in) :: at(3, 3), tau(3, natoms), euler(6)
LOGICAL, INTENT(in) :: lrotation
INTEGER :: i, j
CHARACTER*2 :: label(103)
DATA label /" H", "He", "Li", "Be", " B", " C", " N", " O", " F", "Ne", &
"Na", "Mg", "Al", "Si", " P", " S", "Cl", "Ar", " K", "Ca", &
"Sc", "Ti", " V", "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn", &
"Ga", "Ge", "As", "Se", "Br", "Kr", "Rb", "Sr", " Y", "Zr", &
"Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag", "Cd", "In", "St", &
"Sb", "Te", " I", "Xe", "Cs", "Ba", "La", "Ce", "Pr", "Nd", &
"Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb", &
"Lu", "Hf", "Ta", " W", "Re", "Os", "Ir", "Pt", "Au", "Hg", &
"Tl", "Pb", "Bi", "Po", "At", "Rn", "Fr", "Ra", "Ac", "Th", &
"Pa", " U", "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm", &
"Md", "No", "Lr"/
WRITE(bunit,'("HEADER PROTEIN")')
WRITE(bunit,'("COMPND UNNAMED")')
WRITE(bunit,'("AUTHOR GENERATED BY ...")')
IF (lrotation) &
WRITE(bunit,'("CRYST1",3F9.3,3F7.2,1X,A10,I3)') euler, "P 1", 1
DO i = 1, natoms
WRITE(bunit,'("ATOM ",I5,1X,A2,3X,2A3,I3,3X,F9.3,2F8.3,2F6.2," ")') &
i, label(ityp(i)), "UKN", "", 1, (tau(j, i), j=1,3), 1.0d0, 0.0d0
END DO
WRITE(bunit,'("MASTER 0 0 0 0 0 0 0 0 ", I4," 0 ",I4," 0")') natoms, natoms
WRITE(bunit,'("END")')
RETURN
END SUBROUTINE write_pdb
! PDB File Format
!---------------------------------------------------------------------------
!Field | Column | FORTRAN |
! No. | range | format | Description
!---------------------------------------------------------------------------
! 1. | 1 - 6 | A6 | Record ID (eg ATOM, HETATM)
! 2. | 7 - 11 | I5 | Atom serial number
! - | 12 - 12 | 1X | Blank
! 3. | 13 - 16 | A4 | Atom name (eg " CA " , " ND1")
! 4. | 17 - 17 | A1 | Alternative location code (if any)
! 5. | 18 - 20 | A3 | Standard 3-letter amino acid code for residue
! - | 21 - 21 | 1X | Blank
! 6. | 22 - 22 | A1 | Chain identifier code
! 7. | 23 - 26 | I4 | Residue sequence number
! 8. | 27 - 27 | A1 | Insertion code (if any)
! - | 28 - 30 | 3X | Blank
! 9. | 31 - 38 | F8.3 | Atom's x-coordinate
! 10. | 39 - 46 | F8.3 | Atom's y-coordinate
! 11. | 47 - 54 | F8.3 | Atom's z-coordinate
! 12. | 55 - 60 | F6.2 | Occupancy value for atom
! 13. | 61 - 66 | F6.2 | B-value (thermal factor)
! - | 67 - 67 | 1X | Blank
! 14. | 68 - 68 | I3 | Footnote number
!---------------------------------------------------------------------------
Reference in New Issue View Git Blame Copy Permalink