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

A new code (cpvib.x) for calculating normal modes, Born charge tensors 

and infrared cross-section for isolated molecules/cluster in vacuum.
It uses the CP as an underlying DFT engine, and a the frozen-phonon
approach for calculating the energy Hessian. [silviu]


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@2512 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
silviu 2005-11-24 21:17:29 +00:00
parent 06380219a4
commit cf3931813d
10 changed files with 2904 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

215
VIB/Makefile Normal file
View File

@ -0,0 +1,215 @@
# Makefile for CP/FPMD
include ../make.sys
CPVIB_OBJ = \
vibrations.o \
environment.o \
rdiagh.o \
rdiaghg.o \
read_input_vib.o \
viblib.o
CP_OBJS = \
../CPV/atoms_type.o \
../CPV/band_type.o \
../CPV/bessel.o \
../CPV/berryion.o \
../CPV/bforceion.o \
../CPV/blacs.o \
../CPV/brillouin.o \
../CPV/cg.o \
../CPV/cg_sub.o \
../CPV/cglib.o \
../CPV/chargedensity.o \
../CPV/chargemix.o \
../CPV/charge_types.o \
../CPV/chi2.o \
../CPV/compute_fes_grads.o \
../CPV/compute_scf.o \
../CPV/cp_emass.o \
../CPV/cp_fpmd.o \
../CPV/cp_restart.o \
../CPV/cplib.o \
../CPV/cpr_mod.o \
../CPV/cpr.o \
../CPV/cprsub.o \
../CPV/crayfft.o \
../CPV/dealloc.o \
../CPV/diis.o \
../CPV/dipol_matrix.o \
../CPV/dforceb.o \
../CPV/efermi.o \
../CPV/efield.o \
../CPV/eigsp.o \
../CPV/eigs0.o \
../CPV/electrons.o \
../CPV/emptystates.o \
../CPV/ensemble_dft.o \
../CPV/exch_corr.o \
../CPV/fftdrv.o \
../CPV/fft.o \
../CPV/fields_type.o \
../CPV/fnl.o \
../CPV/forces.o \
../CPV/fromscra.o \
../CPV/greenf.o \
../CPV/grid.o \
../CPV/gsmesh.o \
../CPV/gtable.o \
../CPV/guess.o \
../CPV/init.o \
../CPV/init_run.o \
../CPV/input.o \
../CPV/interfaces_main.o \
../CPV/interfaces.o \
../CPV/ions.o \
../CPV/ions_positions.o \
../CPV/ksstates.o \
../CPV/macdep.o \
../CPV/main.o \
../CPV/mainvar.o \
../CPV/main_loops.o \
../CPV/cplib_meta.o \
../CPV/metaxc.o \
../CPV/modules.o \
../CPV/move_electrons.o \
../CPV/nl_base.o \
../CPV/nlcc.o \
../CPV/nl.o \
../CPV/optical.o \
../CPV/ortho_base.o \
../CPV/ortho.o \
../CPV/para.o \
../CPV/path_routines.o \
../CPV/periodic.o \
../CPV/phasefactor.o \
../CPV/polarization.o \
../CPV/potentials.o \
../CPV/print_out.o \
../CPV/problem_size.o \
../CPV/pseudo_base.o \
../CPV/pseudopot.o \
../CPV/qmatrixd.o \
../CPV/qqberry.o \
../CPV/read_pseudo.o \
../CPV/redis.o \
../CPV/restart.o \
../CPV/restart_sub.o \
../CPV/rsmesh.o \
../CPV/runcg_ion.o \
../CPV/runcg.o \
../CPV/runcp.o \
../CPV/rundiis.o \
../CPV/runsd.o \
../CPV/scalapack.o \
../CPV/smcp.o \
../CPV/smd_modules.o \
../CPV/smd.o \
../CPV/smlam.o \
../CPV/spharmonic.o \
../CPV/spline.o \
../CPV/stop_run.o \
../CPV/stress.o \
../CPV/turbo.o \
../CPV/util.o \
../CPV/vanderwaals.o \
../CPV/cp_version.o \
../CPV/wannier_base.o \
../CPV/wannier.o \
../CPV/waveinit.o \
../CPV/wave.o \
../CPV/wave_types.o \
../CPV/wf.o \
../CPV/$(WRAPPERS)
LOBJS = \
../CPV/adjef.o \
../CPV/entropy.o \
../CPV/forceconv.o \
../CPV/geninv.o \
../CPV/indices.o \
../CPV/miller.o
MODULES = \
../Modules/atom.o \
../Modules/autopilot.o \
../Modules/basic_algebra_routines.o \
../Modules/berry_phase.o \
../Modules/cell_base.o \
../Modules/check_stop.o \
../Modules/clocks.o \
../Modules/coarsegrained_vars.o \
../Modules/constants.o \
../Modules/constraints_module.o \
../Modules/control_flags.o \
../Modules/descriptors.o \
../Modules/electrons_base.o \
../Modules/energies.o \
../Modules/fft_base.o \
../Modules/fft_scalar.o \
../Modules/fft_types.o \
../Modules/functionals.o \
../Modules/griddim.o \
../Modules/input_parameters.o \
../Modules/io_base.o \
../Modules/io_files.o \
../Modules/io_global.o \
../Modules/ions_base.o \
../Modules/ions_nose.o \
../Modules/kind.o \
../Modules/mp_buffers.o \
../Modules/mp_global.o \
../Modules/mp_wave.o \
../Modules/mp.o \
../Modules/parallel_types.o \
../Modules/path_base.o \
../Modules/path_formats.o \
../Modules/path_variables.o \
../Modules/path_opt_routines.o \
../Modules/path_io_routines.o \
../Modules/path_reparametrisation.o \
../Modules/parallel_include.o \
../Modules/parameters.o \
../Modules/parser.o \
../Modules/printout_base.o \
../Modules/pseudo_types.o \
../Modules/ptoolkit.o \
../Modules/read_cards.o \
../Modules/read_namelists.o \
../Modules/readpseudo.o \
../Modules/recvec.o \
../Modules/shmem_include.o \
../Modules/sic.o \
../Modules/smallbox.o \
../Modules/splinelib.o \
../Modules/stick_base.o \
../Modules/uspp.o \
../Modules/version.o \
../Modules/wavefunctions.o \
../Modules/wave_base.o \
../Modules/timestep.o \
../Modules/xml_io_base.o
WRAPPERS = wrapper.o
all : cpvib.x
cpvib: cpvib.x
cpvib.x : $(CP_OBJS) $(CPVIB_OBJ) $(LOBJS) $(LIBOBJS) ../CPV/cp.x vibstart.o
$(MPIF90) $(LDFLAGS) -g -o cpvib.x -I../CPV \
vibstart.o $(CP_OBJS) $(CPVIB_OBJ) $(LOBJS) $(MODULES) \
../flib/eispack.o $(LIBOBJS) $(LIBS)
- ( cd ../bin ; ln -fs ../CPVIB/cpvib.x . )
cpvib_version.o : cpvibver.h
cpvibver.h :
echo "CHARACTER(LEN=70), PARAMETER :: version_date = '"`date`"'" \
> cpvibver.h
clean :
- /bin/rm -f cpvib.x *.o *.mod version.h *.i core* *.F90 fort* \
*.cpp *.d work.pc *.s
include make.depend

212
VIB/environment.f90 Normal file
View File

@ -0,0 +1,212 @@
!
! Copyright (C) 2002-2005 FPMD-CPV groups
! 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 .
!
!==-----------------------------------------------------------------------==!
MODULE environment
!==-----------------------------------------------------------------------==!
USE kinds
USE io_files, ONLY: crash_file, crashunit, &
stop_file, stopunit
IMPLICIT NONE
SAVE
PRIVATE
REAL(DP) :: start_seconds
REAL(DP) :: start_cclock_val
PUBLIC :: environment_start
PUBLIC :: environment_end
PUBLIC :: opening_date_and_time, closing_date_and_time
PUBLIC :: start_cclock_val
!==-----------------------------------------------------------------------==!
CONTAINS
!==-----------------------------------------------------------------------==!
SUBROUTINE environment_start( )
USE io_global, ONLY: stdout, ionode
USE mp_global, ONLY: mpime, nproc
USE parser, ONLY: int_to_char
use para_mod, ONLY: me, node
use mp, only: mp_env
USE cp_version
LOGICAL :: texst
REAL(DP) :: elapsed_seconds, cclock
EXTERNAL elapsed_seconds, cclock
INTEGER :: nchar
CHARACTER(LEN=80) :: uname
CHARACTER(LEN=80) :: version_str
CALL init_clocks( .TRUE. )
CALL start_clock( 'CP' )
start_seconds = elapsed_seconds()
start_cclock_val = cclock( )
version_str = TRIM (version_number) // " - " // TRIM (version_date)
! ... Temporary for para_mod
me = mpime + 1
! ... search for file CRASH and delete it
IF( ionode ) THEN
INQUIRE( FILE=TRIM(crash_file), EXIST=texst )
IF( texst ) THEN
OPEN( UNIT=crashunit, FILE=TRIM(crash_file), STATUS='OLD' )
CLOSE( UNIT=crashunit, STATUS='DELETE' )
END IF
END IF
! ... each processor other than me=1 opens its own standard output file,
! ... this is mainly for debugging, usually only ionode writes to stdout
IF( .NOT. ionode ) THEN
uname = 'out.' // int_to_char( mpime )
nchar = INDEX(uname,' ') - 1
!
! useful for debugging purposes:
! open( unit = stdout, file = uname(1:nchar),status='unknown')
open( unit = stdout, file='/dev/null', status='unknown' )
END IF
! ... Temporary for para_mod
if (me < 10) then
write(node,'(i1,2x)') me
else if (me < 100) then
write(node,'(i2,1x)') me
else if (me < 1000) then
write(node,'(i3)') me
else
call errore('startup','wow, >1000 nodes !!',nproc)
end if
CALL opening_date_and_time( version_str )
#if defined __MPI
WRITE( stdout,100) nproc, mpime
100 FORMAT(3X,'MPI Parallel Build',/,3X,'Tasks =',I5,' This task id =',I5)
#else
WRITE( stdout,100)
100 FORMAT(3X,'Serial Build')
#endif
RETURN
END SUBROUTINE environment_start
!==-----------------------------------------------------------------------==!
SUBROUTINE environment_end( )
USE io_global, ONLY: stdout, ionode
REAL(DP) :: total_seconds
REAL(DP) :: elapsed_seconds
EXTERNAL elapsed_seconds
IF(ionode) THEN
WRITE( stdout,*)
END IF
CALL print_clock( 'CP' )
CALL stop_clock( 'CP' )
CALL closing_date_and_time( )
total_seconds = elapsed_seconds() - start_seconds
IF(ionode) THEN
WRITE( stdout,'(A,F7.1)') ' ELAPSED SECONDS: ', total_seconds
WRITE( stdout,'(A)') ' JOB DONE.'
WRITE( stdout,3335)
END IF
3335 FORMAT('=',78('-'),'=')
RETURN
END SUBROUTINE environment_end
!==-----------------------------------------------------------------------==!
SUBROUTINE opening_date_and_time( version_str )
USE io_global, ONLY: stdout, ionode
CHARACTER(LEN=*), INTENT(IN) :: version_str
CHARACTER(LEN=9) :: cdate, ctime
CHARACTER(LEN=80) :: time_str
CALL date_and_tim(cdate, ctime)
time_str = 'This run was started on: ' // ctime // ' ' // cdate
! ... write program heading
IF(ionode) THEN
WRITE( stdout,3331)
WRITE( stdout,3332) version_str
WRITE( stdout,3331)
WRITE( stdout,3334) time_str
END IF
3331 FORMAT('=',78('-'),'=')
3332 FORMAT( /, 5X,'Ab-initio DFT calculation of normal modes for',/&
& ,5X,'molecules and clusters in vacuum. The CP code is used as',/&
& ,5x,'the undelying DFT egine.'//&
& ,5X,'Version: ',A60,/&
& ,5x,'Authors of the normal-modes code: Silviu Zilberman',//&
& ,5X,'Authors of the CP code: Alfredo Pasquarello, Kari Laasonen, Andrea Trave, Roberto Car,',/&
& ,5X,' Paolo Giannozzi, Nicola Marzari, Carlo Cavazzoni, Guido Chiarotti,',/&
& ,5X,' Sandro Scandolo, Paolo Focher, Gerardo Ballabio, and others',/)
3334 FORMAT(/,3X,A60,/)
RETURN
END SUBROUTINE opening_date_and_time
!==-----------------------------------------------------------------------==!
SUBROUTINE closing_date_and_time( )
USE io_global, ONLY: stdout, ionode
CHARACTER(LEN=9) :: cdate, ctime
CHARACTER(LEN=80) :: time_str
CALL date_and_tim(cdate, ctime)
time_str = 'This run was terminated on: ' // ctime // ' ' // cdate
IF( ionode ) THEN
WRITE( stdout,*)
WRITE( stdout,3334) time_str
WRITE( stdout,3335)
END IF
3334 FORMAT(3X,A60,/)
3335 FORMAT('=',78('-'),'=')
RETURN
END SUBROUTINE closing_date_and_time
!==-----------------------------------------------------------------------==!
END MODULE environment
!==-----------------------------------------------------------------------==!

101
VIB/rdiagh.f90 Normal file
View File

@ -0,0 +1,101 @@
!
! Copyright (C) 2001-2005 PWSCF 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 .
!
!#include "f_defs.h"
!
!----------------------------------------------------------------------------
SUBROUTINE rdiagh( n, h, ldh, e, v )
!----------------------------------------------------------------------------
!
! ... calculates all the eigenvalues and eigenvectors of a complex
! ... hermitean matrix H . On output, the matrix is unchanged
!
!USE kinds, ONLY : 8
!USE mp_global, ONLY : npool, me_pool, root_pool, intra_pool_comm, my_image_id
!USE mp, ONLY : mp_bcast
!
IMPLICIT NONE
!
! ... on INPUT
!
INTEGER :: n, & ! dimension of the matrix to be diagonalized
ldh ! leading dimension of h, as declared in
! the calling pgm unit
REAL (KIND=8) :: h(ldh,n) ! matrix to be diagonalized
!
! ... on OUTPUT
!
REAL (KIND=8) :: e(n) ! eigenvalues
REAL (KIND=8) :: v(ldh,n) ! eigenvectors (column-wise)
!
!
!CALL start_clock( 'rdiagh' )
!
CALL rdiagh_lapack( )
!
!CALL stop_clock( 'rdiagh' )
!
RETURN
!
CONTAINS
!
! ... internal procedures
!
!-----------------------------------------------------------------------
SUBROUTINE rdiagh_lapack( )
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
! ... local variables (LAPACK version)
!
INTEGER :: lwork, nb, info
INTEGER, EXTERNAL :: ILAENV
! ILAENV returns optimal block size "nb"
REAL (KIND=8), ALLOCATABLE :: work(:)
!
!
! ... check for the block size
!
nb = ILAENV( 1, 'DSYTRD', 'U', n, - 1, - 1, - 1 )
!
IF ( nb < 1 ) nb = MAX( 1, n )
!
lwork = ( nb + 3 ) * n
!
! ... only the first processor diagonalize the matrix
!
!IF ( me_pool == root_pool ) THEN
!
! ... allocate workspace
!
v = h
!
ALLOCATE( work( lwork ) )
!
CALL DSYEV( 'V', 'U', n, v, ldh, e, work, lwork, info )
!
if (abs(info).ne.0) then
!CALL errore( 'rdiagh', 'info =/= 0', ABS( info ) )
write (6,*) 'error: rdiagh - info =/= 0 ', ABS( info )
end if
!
! ... deallocate workspace
!
DEALLOCATE( work )
!
!END IF
!
!CALL mp_bcast( e, root_pool, intra_pool_comm )
!CALL mp_bcast( v, root_pool, intra_pool_comm )
!
RETURN
!
END SUBROUTINE rdiagh_lapack
!
END SUBROUTINE rdiagh

153
VIB/rdiaghg.f90 Normal file
View File

@ -0,0 +1,153 @@
!
! Copyright (C) 2005 PWSCF 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 .
!
!#include "f_defs.h"
!
!----------------------------------------------------------------------------
SUBROUTINE rdiaghg( n, m, h, s, ldh, e, v )
!----------------------------------------------------------------------------
!
! ... calculates eigenvalues and eigenvectors of the generalized problem
! ... Hv=eSv, with H symmetric matrix, S overlap matrix .
! ... On output both matrix are unchanged
! ... Uses LAPACK routines
!
!USE kinds, ONLY : 8
!USE mp, ONLY : mp_bcast
!USE mp_global, ONLY : npool, me_pool, root_pool, intra_pool_comm, my_image_id
!
IMPLICIT NONE
!
! ... on INPUT
!
INTEGER :: n, m, ldh
! dimension of the matrix to be diagonalized
! number of eigenstates to be calculated
! leading dimension of h, as declared in the calling pgm unit
REAL(KIND=8) :: h(ldh,n)
! matrix to be diagonalized
REAL(KIND=8) :: s(ldh,n)
! overlap matrix
!
! ... on OUTPUT
!
REAL(KIND=8) :: e(n)
! eigenvalues
REAL(KIND=8) :: v(ldh,m)
! eigenvectors (column-wise)
!
! ... LOCAL variables
!
INTEGER :: lwork, nb, ILAENV, mm, info
! ILAENV returns optimal block size "nb"
! mm = number of calculated eigenvectors
EXTERNAL ILAENV
INTEGER, ALLOCATABLE :: iwork(:), ifail(:)
REAL(KIND=8), ALLOCATABLE :: sdum(:,:), hdum(:,:), work(:)
LOGICAL :: all_eigenvalues
!
!
!CALL start_clock( 'cdiaghg' )
!
all_eigenvalues = ( m == n )
!
! ... check for optimal block size
!
nb = ILAENV( 1, 'DSYTRD', 'U', n, -1, -1, -1 )
!
IF ( nb < 1 .OR. nb >= n ) THEN
!
lwork = 8 * n
!
ELSE
!
lwork = ( nb + 3 ) * n
!
END IF
!
! ... allocate workspace
!
ALLOCATE( work( lwork ) )
ALLOCATE( sdum( ldh, n ) )
!
IF ( .NOT. all_eigenvalues ) THEN
!
ALLOCATE( hdum( ldh, n ) )
ALLOCATE( iwork( 5*n ) )
ALLOCATE( ifail( n ) )
!
END IF
!
! ... input s and (see below) h are copied so that they are not destroyed
!
sdum = s
!
! ... only the first processor diagonalize the matrix
!
!IF ( me_pool == root_pool ) THEN
!
IF ( all_eigenvalues ) THEN
!
! ... calculate all eigenvalues
!
v(:,1:n) = h(:,:)
!
#if defined (__AIX)
!
! ... there is a name conflict between essl and lapack ...
!
CALL DSYGV( 1, v, ldh, sdum, ldh, e, v, ldh, n, work, lwork )
!
info = 0
!
#else
CALL DSYGV( 1, 'V', 'U', n, v, ldh, sdum, ldh, e, work, &
lwork, info )
!
#endif
ELSE
!
! ... calculate only m lowest eigenvalues
!
hdum = h
!
CALL DSYGVX( 1, 'V', 'I', 'U', n, hdum, ldh, sdum, ldh, &
0.D0, 0.D0, 1, m, 0.D0, mm, e, v, ldh, work, lwork, &
iwork, ifail, info )
!
END IF
!
!CALL errore( 'rdiaghg', 'info =/= 0', ABS( info ) )
if (info.ne.0) then
write (6,*) 'error: rdiaghg - info =/= 0i ',ABS( info )
end if
!
!END IF
!
! ... broadcast eigenvectors and eigenvalues to all other processors
!
!CALL mp_bcast( e, root_pool, intra_pool_comm )
!CALL mp_bcast( v, root_pool, intra_pool_comm )
!
! ... deallocate workspace
!
IF ( .NOT. all_eigenvalues ) THEN
!
DEALLOCATE( ifail )
DEALLOCATE( iwork )
DEALLOCATE( hdum )
!
END IF
!
DEALLOCATE( sdum )
DEALLOCATE( work )
!
!CALL stop_clock( 'cdiaghg' )
!
RETURN
!
END SUBROUTINE rdiaghg

138
VIB/read_input_vib.f90 Normal file
View File

@ -0,0 +1,138 @@
!
! Copyright (C) 2001-2005 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 .
!
#include "f_defs.h"
!
!----------------------------------------------------------------------------
SUBROUTINE read_input_vib()
!----------------------------------------------------------------------------
!
! This routine reads the control variables for the program cpvib.
! from prefix.vib.inp input file.
!
! ... modules
!
USE io_files, ONLY : outdir, prefix
USE io_global, ONLY : ionode, ionode_id, stdout
USE kinds, ONLY : DP
USE mp, ONLY : mp_bcast
USE vibrations, ONLY : delta, save_freq, &
trans_inv_max_iter, trans_inv_conv_thr, vib_restart_mode, &
trans_inv_flag, trans_rot_inv_flag, animate
!
IMPLICIT NONE
!
! ... local variables
!
INTEGER :: ios, dirlen
CHARACTER (LEN=256) :: restart_label, input_filename
LOGICAL :: inp_file_exists
!
NAMELIST / INPUT_VIB / delta, save_freq, trans_inv_max_iter, &
trans_inv_conv_thr, restart_label, trans_inv_flag,trans_rot_inv_flag, &
animate
!
IF ( ionode ) THEN
!
! Set input file name
!
dirlen = INDEX(outdir,' ') - 1
input_filename = TRIM(prefix)//'.vib.inp'
input_filename = outdir(1:dirlen) // '/' // input_filename
WRITE (stdout,*) 'Input information is read from:', input_filename
!
! ... Open input file
!
INQUIRE (FILE = input_filename, EXIST = inp_file_exists)
IF ( .NOT. inp_file_exists) &
CALL errore( 'read_input_vib', 'input file absent ', 1 )
OPEN( 99 , FILE = input_filename , STATUS = 'old' )
!
! ... set default values for variables in namelist
!
delta = 0.05
save_freq = 1
trans_inv_flag = .TRUE.
trans_rot_inv_flag = .TRUE.
trans_inv_max_iter = 50
trans_inv_conv_thr = 1.0D-15
vib_restart_mode = 2
restart_label = 'auto'
animate = .FALSE.
!
! ... reading the namelist input_vib
!
READ( 99, input_vib, ERR = 200, IOSTAT = ios )
!
200 CALL errore( 'readin_input_vib', 'reading input_vib namelist', ABS(ios) )
!
! ... Check all namelist variables
!
IF (delta <= 0.D0) &
CALL errore (' readin_input_vib', ' Wrong delta ', 1)
!
IF (trans_inv_max_iter <= 0.D0) &
CALL errore (' readin_input_vib', ' Wrong trans_inv_max_iter ', 1)
!
IF (trans_inv_conv_thr <= 0.D0) &
CALL errore (' readin_input_vib', ' Wrong trans_inv_conv_thr ', 1)
!
SELECT CASE ( TRIM( restart_label ) )
CASE( 'from_scratch' )
vib_restart_mode = 0
CASE ( 'restart' )
vib_restart_mode = 1
CASE ( 'auto' )
vib_restart_mode = 2
CASE DEFAULT
CALL errore( 'read_input_vib ', &
'unknown vib_restart_mode ' // TRIM( restart_label ), 1 )
END SELECT
END IF
!
! ... Broadcasting input variables
!
CALL mp_bcast( delta, ionode_id )
CALL mp_bcast( save_freq, ionode_id )
CALL mp_bcast( trans_inv_max_iter, ionode_id )
CALL mp_bcast( trans_inv_conv_thr, ionode_id )
CALL mp_bcast( trans_inv_flag, ionode_id )
CALL mp_bcast( trans_rot_inv_flag, ionode_id )
CALL mp_bcast( vib_restart_mode, ionode_id )
CALL mp_bcast( animate, ionode_id )
!
! ... Printing input data to stdout
!
IF (ionode) THEN
WRITE (stdout,*) '------------------------------------------------------'
WRITE (stdout,*) ' SUMMERY OF VIBRATIONAL INPUT VARIABLES: '
WRITE (stdout,*) ''
WRITE (stdout,101) delta
WRITE (stdout,102) save_freq
WRITE (stdout,103) trans_inv_flag
WRITE (stdout,104) trans_inv_max_iter
WRITE (stdout,105) trans_inv_conv_thr
WRITE (stdout,106) trans_rot_inv_flag
WRITE (stdout,107) restart_label
WRITE (stdout,108) animate
WRITE (stdout,*) '------------------------------------------------------'
!
101 FORMAT(3x,'delta = ',3x,f10.4)
102 FORMAT(3x,'save_freq = ',3x,I10)
103 FORMAT(3x,'trans_inv_flag = ',3x,L10)
104 FORMAT(3x,'trans_inv_max_iter = ',3x,I10)
105 FORMAT(3x,'trans_inv_conv_thr = ',3x,D10.4)
106 FORMAT(3x,'trans_rot_inv_flag = ',3x,L10)
107 FORMAT(3x,'restart_label = ',3x,A10)
108 FORMAT(3x,'animate = ',3x,L10)
!
CLOSE(99)
END IF
!
RETURN
!
END SUBROUTINE read_input_vib

328
VIB/viblib.f90 Normal file
View File

@ -0,0 +1,328 @@
!
! Copyright (C) 2001-2005 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 .
!
!-------------------------------------------------------------------------
!
! Auxiliary functions for the vibrational analysis
!
!-------------------------------------------------------------------------
SUBROUTINE print_matrix (matrix, dim1, dim2, title, filep)
!
! ... modules
!
USE kinds, ONLY: DP
!
IMPLICIT NONE
!
! ... input variables
!
INTEGER, INTENT(IN) :: dim1, dim2, filep
CHARACTER (len=*), INTENT(IN) :: title
REAL (KIND=DP), INTENT(IN) :: matrix(dim1,dim2)
!
! ... local variables
!
INTEGER :: i,j
!
! ...
!
WRITE (filep,*)
WRITE (filep,*) TRIM(title)
WRITE (filep,*)
DO i=1,dim1
DO j=1,dim2
WRITE (filep,FMT='(F12.5,3X)',ADVANCE='NO') matrix(i,j)
END DO
WRITE (filep,*)
END DO
WRITE (filep,*)
!
RETURN
END SUBROUTINE print_matrix
!
!
! ----------------------------------------------
!
!
SUBROUTINE write_xyz(tau,free_text,position_in_file,file_name)
!
! Printout of the coordinates in an xyz format
! tau - coordinates
! free_text - text for the second line of output
! position_in_file - either 'APPEND' (for an xyz animation),
! 'REWIND' to overwrite, or 'ASIS'
!
! ... modules
!
USE constants, ONLY : bohr_radius_angs
USE input_parameters, ONLY : atom_label
USE io_files, ONLY : iunxyz
USE ions_base, ONLY : nat, ityp
USE kinds, ONLY : DP
!
IMPLICIT NONE
!
! ... input variables
!
CHARACTER (LEN=*), INTENT(IN) :: position_in_file, free_text
CHARACTER (LEN=*), INTENT(IN) :: file_name
REAL (KIND=DP), INTENT(IN) :: tau(3,nat)
!
! ... local variables
!
INTEGER :: atom
CHARACTER (LEN=*), PARAMETER :: xyz_fmt = "(A2,3(2X,F14.10))"
OPEN( UNIT = iunxyz, FILE = file_name, &
STATUS = "UNKNOWN", ACTION = "WRITE", &
POSITION = TRIM(position_in_file))
!
WRITE( UNIT = iunxyz, FMT = '(I5,/)' ) nat
WRITE( UNIT = iunxyz, FMT = '(A ,/)' ) free_text
!
DO atom = 1, nat
!
WRITE( UNIT = iunxyz, FMT = xyz_fmt ) &
TRIM( atom_label( ityp( atom ) ) ), &
tau(1,atom)*bohr_radius_angs, &
tau(2,atom)*bohr_radius_angs, &
tau(3,atom)*bohr_radius_angs
!
END DO
!
CLOSE(iunxyz)
!
RETURN
END SUBROUTINE write_xyz
!
!
! ----------------------------------------------
!
!
SUBROUTINE calculate_dipole (dipole, dipole_moment,tau)
!
! Driver routine for the calculation of the dipole moment.
! Currently it uses only subroutine 'poles' from cplib.f90
! It could easily be extended to use Wannier functions as well.
!
! ... modules
!
USE cp_main_variables, ONLY : irb, eigrb, bec, rhor, rhog, rhos
USE electrons_base, ONLY : nspin
USE energies, ONLY : ekin, enl
USE ions_base, ONLY : nat
USE kinds, ONLY : DP
USE wavefunctions_module, ONLY : c0
USE uspp, ONLY : becsum
!
IMPLICIT NONE
!
! ... input variables
!
REAL (KIND=DP), INTENT(IN) :: tau(3,nat)
!
! ... output variables
!
REAL (KIND=DP), INTENT(OUT) :: dipole(3), dipole_moment
!
! ... local variables
!
REAL (KIND=DP) :: quadrupole
INTEGER :: nfi=10 ! ... dummy value
LOGICAL :: ion_flag = .TRUE., coc_flag=.TRUE.
!
!
!
dipole = 0.0
dipole_moment = 0.0
rhor = 0.0
rhog = 0.0
rhos = 0.0
!
CALL rhoofr(nfi,c0,irb,eigrb,bec,becsum,rhor,rhog,rhos,enl,ekin)
IF (nspin.EQ.2) THEN
rhor(:,1)=rhor(:,1)+rhor(:,2)
END IF
CALL poles (dipole_moment,dipole,quadrupole,rhor(:,1),&
ion_flag,tau,coc_flag)
!
! STILL HAVE TO ADD IN THE WANNIER BASED DIPOLE
! AS AN ALTERNATIVE WAY OF DOING THE CALCULATION
!
RETURN
END SUBROUTINE calculate_dipole
!
!
! ----------------------------------------------
!
!
SUBROUTINE orthonormalize(M,dim1,dim2)
!
! modified Gram-Schimd - orthogonalizing for
! dim1 vectors of length dim2
!
!
! ... modules
!
USE kinds, ONLY : DP
!
IMPLICIT NONE
!
!... input variables
!
INTEGER, INTENT(IN) :: dim1,dim2 !dimensions of input matrix
REAL(KIND=DP),INTENT(INOUT):: M(dim1,dim2) !input matrix
!
!... local variables
!
REAL (KIND=DP) :: t(dim2,dim2), z
INTEGER :: i,j,k
!
t=0.0
!
DO k=1,dim2
!
! ... normalizing vector k
z=0.0
DO i=1,dim1
z=z+M(i,k)**2
END DO
t(k,k)=SQRT(z)
M(:,k)=M(:,k)/t(k,k)
!
! ... removing the projection of vector k
! ... from the remaining vectors k+1...dim2
!
DO j=k+1,dim2
z=0.0
DO i=1,dim1
z=z+M(i,j)*M(i,k)
END DO
!
t(k,j)=z
!
DO i=1,dim1
M(i,j)=M(i,j)-t(k,j)*M(i,k)
END DO
END DO
!
END DO
!
RETURN
END SUBROUTINE orthonormalize
!
!
! ----------------------------------------------
!
!
SUBROUTINE symmetrize_matrix(M, dim)
!
! IMPOSING SYMMETRY ON MATRIX M
! M(i,j)==M(j,i)
!
! ... modules
!
USE kinds, ONLY : DP
!
IMPLICIT NONE
!
! ... input/output variables
!
INTEGER, INTENT(IN) :: dim
REAL (KIND=DP), INTENT(INOUT) :: M(dim,dim)
!
! ... local variables
!
INTEGER :: i,j
REAL (KIND=DP) :: sum
!
!
DO i=1,dim
DO j = i+1,dim
sum = (M(i,j)+M(j,i))/2
M(i,j) = sum
M(j,i) = sum
END DO
END DO
!
RETURN
END SUBROUTINE symmetrize_matrix
!
!
! ----------------------------------------------
!
!
SUBROUTINE relax_wavefunction (fion)
!
! A DRIVER ROUTINE FOR WAVE FUNCTION RELAXATION
!
! ... modules
!
USE cell_base, ONLY : b1, b2, b3
USE cg_module, ONLY : tcg
USE control_flags, ONLY : tprnfor, thdyn
USE cp_electronic_mass, ONLY : emass
USE cp_main_variables, ONLY : eigr, ei1, ei2, ei3, &
sfac, irb, eigrb, taub, nfi
USE efield_module, ONLY : tefield, efield_update
USE electrons_nose, ONLY : fccc
USE energies, ONLY : eself, etot
USE from_scratch_module, ONLY : from_scratch
USE gvecs, ONLY : ngs
USE ions_positions, ONLY : tau0
USE kinds, ONLY : DP
USE parameters, ONLY : natx
USE phase_factors_module, ONLY : strucf
USE reciprocal_vectors, ONLY : mill_l
USE time_step, ONLY : dt2
!
IMPLICIT NONE
!
! ... input variables
!
REAL (KIND=DP) :: fion(3,natx)
!
! ... local variables
!
LOGICAL :: tfirst, tlast
REAL (KIND=DP) :: enthal
REAL (KIND=DP) :: dt2bye, enb, enbi, ccc
!
! ... for smooth restart in the new coordinates
!
dt2bye = dt2 / emass
CALL initbox( tau0, taub, irb )
CALL phbox( taub, eigrb )
CALL phfac( tau0, ei1, ei2, ei3, eigr )
CALL strucf( sfac, ei1, ei2, ei3, mill_l, ngs )
IF ( thdyn ) CALL formf( tfirst, eself )
IF (tefield ) CALL efield_update( eigr )
!
!... relax wavefunction in new position
!
IF (tcg) THEN
tprnfor = .TRUE. ! ... atomic forces are calculated only at the
! end of the wavefunction relaxation process
CALL move_electrons( nfi, tfirst, tlast, b1, b2, b3, fion, &
enthal, enb, enbi, fccc, ccc, dt2bye )
tprnfor = .FALSE.
ELSE
tprnfor = .FALSE. ! ... no need to calculate atomic forces at
! each step of the MD damped dynamics
CALL cprmain ( tau0, fion, etot )
!
! ... one more scf step to get the forces on the nuclei
!
tprnfor = .TRUE.
CALL move_electrons( nfi, tfirst, tlast, b1, b2, b3, fion, &
enthal, enb, enbi, fccc, ccc, dt2bye )
tprnfor = .FALSE.
END IF
!
RETURN
END SUBROUTINE relax_wavefunction

1540
VIB/vibrations.f90 Normal file
View File

File diff suppressed because it is too large Load Diff

213
VIB/vibstart.f90 Normal file
View File

@ -0,0 +1,213 @@
!
! Copyright (C) 2002-2005 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 .
!
#include "f_defs.h"
!
!==============================================================================
! Program CPVIB (programmed by Silviu Zilberman)
!
! Calculates the vibrational modes, Born effective charges and infrared
! cross-section for isolated molecules/clusters in vacuum.
! The programs uses the CP code as the underlying DFT engine
! for wave function relaxation (SCF) and the frozen phonon method for
! construction of the energy hessian (dynamical matrix),
! i.e. displace each atom in all three cartesian directions to construct
! the numerical first derivative of the forces.
!
! The program uses two mandatory input files (optional files are discussed
! further below):
! 1. The regular CP input file, which presumably was used to generate a
! ground state nuclear configuration (local minimum).
! 2. A file "prefix.vib.inp", where prefix is defined in the CP input
! file.
!
! prefix.vib.inp contains one name list:
!
! &INPUT_VIB
! ...
! /
!
! Parameters:
!
! delta real ( default = 0.05 )
! displacement step, a uniform value for all atoms
!
! save_freq integer ( default = 1 )
! Save restart information every save_freq displacements.
! Since it is usually a small file, containing mainly
! the intermediate hessian, it is suggested not to
! change the default.
! restart_label character ( default = 'auto' )
! 'from_scratch': starts a new calculation
! 'restart' : continue a previous calculation
! a file 'prefix.vib.restart' must be present
! 'auto' : if a file 'prefix.vib.restart' is present
! the program will pick it up and continue
! the previous calculation, or if absent,
! it will start a new calculation
!
! trans_inv_flag logical ( default = .true. )
! impose translational invariance on the energy hessian
! to ensure that the modes associated with a rigid
! translation would have zero frequency
!
! trans_inv_max_iter integer ( default = 50 )
! the maximal number of steps in the iterative procedure
! that removes the rigid translation modes
!
! trans_inv_conv_thr real ( default = 1.0D-15 )
! the threshhold for convergence of the iterative procedure
! for removal of the rigid translation modes
!
! trans_rot_inv_flag logical ( default = .true. )
! remove also the rigid rotation modes
!
! animate logical ( default = .false. )
! generate xyz animation files for all normal modes
!
! Additional information about the input/output files:
!
! 1. The program will look for the presence of 'prefix.vib.isotope' file.
! It is a text file containing NAT lines, where NAT is the number of atoms.
! Each line contains a single number, the mass of this particular atom
! in the same order as in the CP input file. The masses are in AMU, i.e.
! the hydrogen is ~1, carbon is 12 etc. If this file is present, these masses
! are used in the calculation. If not, the default masses are used,
! as specified in the CP input, and an isotopes file is created.
! This is useful for testing the isotope shifts, without having to
! redo the full calculation.
!
! 2. After the energy hessian is calculated, the program will write
! the analysis results to 'prefix.vib.analysis'.
! The minimal output is:
! * raw (but symmetrized) hessian
! * diagonal mass matrix
! * harmonic frequencies and the associated effective
! spring constant and effective mode mass
! * raw Born charge tensors for each atom
! * the sum over the Born charges (minus total system charge)
! this quantity should ideally be zero (acoustic sum rule).
! * infrared intensity of each normal mode.
!
! If trans_inv_flag==.true. then the same information as above is repeated
! but this time with imposed translational invariance on the energy
! hessian. In addition, the Born charges are symmetrized and
! the acoustic sum rule is imposed.
!
! if trans_rot_inv_flag==.true. then also rotational modes are projected
! out.
!
! A brief description of the methods:
!
! 1. The energy hessian is constructed by a simple difference formula for
! the numerical first derivative of the forces.
!
! 2. The condition for translational invariance of the hessian is described
! (for instance) in Gonze and Lee, Phys. Rev. B 55(16), 1997, 10355-10368.
! However, simple enforcement of translational invariance breaks the
! symmetry of the hessian. Here symmetrization and translational
! invariance are repeatedly applied until convergence, i.e., until
! the largest change in any matrix element of the hessian is smaller
! than trans_inv_conv_thr parameter.
! A more sophisticated and general algorithm was suggested by
! Nicolas Mounet and Nicola Marzari (Ref. ???) and is already implemented
! in QE.
!
! 3. The algorithm for removal of rotational rigid modes follows closely the
! one implemented in Gaussian03 electronic structure package, as described
! in their white papers (www.gaussian.com).
!
!
! ==========================================================================
! ==== units and constants ====
! ==== ====
! ==== 1 hartree = 1 a.u. ====
! ==== 1 bohr radius = 1 a.u. = 0.529167 Angstrom ====
! ==== 1 rydberg = 1/2 a.u. ====
! ==== 1 electron volt = 1/27.212 a.u. ====
! ==== 1 kelvin *k-boltzm. = 1/(27.212*11606) a.u.'='3.2e-6 a.u ====
! ==== 1 second = 1/(2.4189)*1.e+17 a.u. ====
! ==== 1 proton mass = 1822.89 a.u. ====
! ==== 1 tera = 1.e+12 ====
! ==== 1 pico = 1.e-12 ====
! ==== 1 Volt / meter = 1/(5.1412*1.e+11) a.u. ====
! ==========================================================================
!
!----------------------------------------------------------------------------
PROGRAM cpvib
!----------------------------------------------------------------------------
!
USE control_flags, ONLY : program_name
USE environment, ONLY : environment_start
USE input, ONLY : read_input_file, iosys_pseudo, iosys
USE io_global, ONLY : io_global_start, io_global_getionode
USE kinds, ONLY : DP
USE mp_global, ONLY : mp_global_start
USE mp, ONLY : mp_end, mp_start, mp_env
USE vibrations, ONLY : start_vibrations, calc_hessian, &
analysis, end_vibrations
!
IMPLICIT NONE
!
INTEGER :: mpime, nproc, gid, ionode_id, &
restart_cyc_counter
INTEGER, PARAMETER :: root = 0
LOGICAL :: ionode
REAL (KIND=DP) :: E_minus, dip_minus(3)
!
! ... program starts here
!
program_name = 'CP90'
!
! ... initialize MPI (parallel processing handling)
!
CALL mp_start()
CALL mp_env( nproc, mpime, gid )
CALL mp_global_start( root, mpime, gid, nproc )
!
! ... mpime = processor number, starting from 0
! ... nproc = number of processors
! ... gid = group index
! ... root = index of the root processor
!
! ... initialize input output
!
CALL io_global_start( mpime, root )
CALL io_global_getionode( ionode, ionode_id )
!
CALL environment_start( )
!
! ... readin the input file
!
CALL read_input_file()
!
! ... copy pseudopotential input parameter into internal variables
! ... and read in pseudopotentials and wavefunctions files
!
CALL iosys_pseudo()
!
! ... copy-in input parameters from input_parameter module
!
CALL iosys()
!
CALL init_run()
!
! CALL memstat( 1 )
!
CALL read_input_vib()
CALL start_vibrations(restart_cyc_counter, E_minus, dip_minus)
CALL calc_hessian (restart_cyc_counter, E_minus, dip_minus)
CALL analysis ()
CALL end_vibrations()
!
CALL terminate_run()
!
CALL stop_run( .TRUE. )
!
STOP
!
END PROGRAM cpvib

3
install/makedeps.sh
View File

@ -6,7 +6,7 @@ cd `echo $0 | sed 's/\(.*\)\/.*/\1/'` # extract pathname
TOPDIR=`pwd`
for DIR in Modules clib PW CPV flib pwtools upftools PP PWCOND \
Gamma PH D3 atomic Nmr
Gamma PH D3 atomic Nmr CPVIB
do
# set inter-directory dependencies
case $DIR in
@ -17,6 +17,7 @@ do
PP | PWCOND | Gamma | PH )
DEPENDS="../include ../Modules ../PW ../iotk/src" ;;
D3 | Nmr) DEPENDS="../include ../Modules ../PW ../PH ../iotk/src" ;;
CPVIB ) DEPENDS="../include ../Modules ../iotk/src ../CPV"
esac
# generate dependencies file

3
makedeps.sh
View File

@ -6,7 +6,7 @@ cd `echo $0 | sed 's/\(.*\)\/.*/\1/'` # extract pathname
TOPDIR=`pwd`
for DIR in Modules clib PW CPV flib pwtools upftools PP PWCOND \
Gamma PH D3 atomic Nmr
Gamma PH D3 atomic Nmr CPVIB
do
# set inter-directory dependencies
case $DIR in
@ -17,6 +17,7 @@ do
PP | PWCOND | Gamma | PH )
DEPENDS="../include ../Modules ../PW ../iotk/src" ;;
D3 | Nmr) DEPENDS="../include ../Modules ../PW ../PH ../iotk/src" ;;
CPVIB ) DEPENDS="../include ../Modules ../iotk/src ../CPV"
esac
# generate dependencies file