2005-11-25 05:17:29 +08:00
|
|
|
!
|
|
|
|
! 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
|
2005-12-02 06:40:24 +08:00
|
|
|
WRITE (filep,FMT='(F12.6,3X)',ADVANCE='NO') matrix(i,j)
|
2005-11-25 05:17:29 +08:00
|
|
|
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))
|
|
|
|
!
|
2005-11-26 07:39:31 +08:00
|
|
|
WRITE( UNIT = iunxyz, FMT = '(I5)' ) nat
|
|
|
|
WRITE( UNIT = iunxyz, FMT = '(A)' ) free_text
|
2005-11-25 05:17:29 +08:00
|
|
|
!
|
|
|
|
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
|
|
|
|
!
|
2005-12-29 16:38:27 +08:00
|
|
|
USE kinds, ONLY : DP
|
|
|
|
USE ions_base, ONLY : nat
|
2005-12-19 21:31:56 +08:00
|
|
|
#ifdef DFT_CP
|
2005-12-29 16:38:27 +08:00
|
|
|
USE cp_main_variables, ONLY : irb, eigrb, bec, charge_density=>rhor, &
|
|
|
|
rhog, rhos
|
2005-11-25 05:17:29 +08:00
|
|
|
USE electrons_base, ONLY : nspin
|
|
|
|
USE energies, ONLY : ekin, enl
|
|
|
|
USE wavefunctions_module, ONLY : c0
|
|
|
|
USE uspp, ONLY : becsum
|
2005-12-29 16:38:27 +08:00
|
|
|
USE grid_dimensions, ONLY : nnrx
|
|
|
|
#endif
|
|
|
|
#ifdef DFT_PW
|
|
|
|
USE scf, ONLY : charge_density=>rho
|
|
|
|
USE lsda_mod, ONLY : nspin
|
|
|
|
USE gvect, ONLY : nrxx
|
|
|
|
USE cell_base, ONLY : alat
|
|
|
|
#endif
|
2005-11-25 05:17:29 +08:00
|
|
|
!
|
|
|
|
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.
|
2005-12-29 16:38:27 +08:00
|
|
|
#ifdef DFT_PW
|
|
|
|
REAL(DP) :: tmp_rho(nrxx,nspin)
|
|
|
|
#endif
|
|
|
|
#ifdef DFT_CP
|
|
|
|
REAL(DP) :: tmp_rho(nnrx,nspin)
|
|
|
|
#endif
|
|
|
|
!
|
|
|
|
!
|
|
|
|
!
|
|
|
|
! ... initiating variables
|
|
|
|
dipole = 0.0
|
|
|
|
dipole_moment = 0.0
|
|
|
|
!
|
|
|
|
#ifdef DFT_CP
|
|
|
|
!
|
|
|
|
! ... CP charge density
|
|
|
|
!
|
|
|
|
rhog = 0.0
|
|
|
|
rhos = 0.0
|
2005-11-25 05:17:29 +08:00
|
|
|
!
|
2005-12-29 16:38:27 +08:00
|
|
|
CALL rhoofr(nfi,c0,irb,eigrb,bec,becsum,charge_density,rhog,rhos,enl,ekin)
|
|
|
|
#endif
|
|
|
|
|
|
|
|
#ifdef DFT_PW
|
2005-11-25 05:17:29 +08:00
|
|
|
!
|
2005-12-29 16:38:27 +08:00
|
|
|
! ... PW charge density
|
2005-11-25 05:17:29 +08:00
|
|
|
!
|
2005-12-29 16:38:27 +08:00
|
|
|
#endif
|
2005-11-25 05:17:29 +08:00
|
|
|
!
|
2005-12-29 16:38:27 +08:00
|
|
|
tmp_rho = charge_density
|
2005-11-25 05:17:29 +08:00
|
|
|
IF (nspin.EQ.2) THEN
|
2005-12-29 16:38:27 +08:00
|
|
|
tmp_rho(:,1)=tmp_rho(:,1)+tmp_rho(:,2)
|
2005-11-25 05:17:29 +08:00
|
|
|
END IF
|
2005-12-29 16:38:27 +08:00
|
|
|
!
|
|
|
|
! ... computing dipole
|
|
|
|
!
|
|
|
|
#ifdef DFT_CP
|
|
|
|
CALL poles (dipole_moment,dipole,quadrupole,tmp_rho(:,1),&
|
2005-11-25 05:17:29 +08:00
|
|
|
ion_flag,tau,coc_flag)
|
2005-12-29 16:38:27 +08:00
|
|
|
#endif
|
|
|
|
#ifdef DFT_PW
|
|
|
|
CALL poles (dipole_moment,dipole,quadrupole,tmp_rho(:,1),&
|
|
|
|
ion_flag,tau*alat,coc_flag)
|
|
|
|
#endif
|
2005-11-25 05:17:29 +08:00
|
|
|
!
|
|
|
|
! STILL HAVE TO ADD IN THE WANNIER BASED DIPOLE
|
2005-12-29 16:38:27 +08:00
|
|
|
! AS AN ALTERNATIVE WAY FOR CALCULATING THE DIPOLE
|
2005-11-25 05:17:29 +08:00
|
|
|
!
|
|
|
|
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
|
2005-12-23 07:35:20 +08:00
|
|
|
!
|
|
|
|
!
|
|
|
|
! ----------------------------------------------
|
|
|
|
!
|
|
|
|
!
|
2005-11-25 05:17:29 +08:00
|
|
|
SUBROUTINE relax_wavefunction (fion)
|
|
|
|
!
|
|
|
|
! A DRIVER ROUTINE FOR WAVE FUNCTION RELAXATION
|
|
|
|
!
|
2005-12-19 21:31:56 +08:00
|
|
|
!
|
2005-11-25 05:17:29 +08:00
|
|
|
! ... modules
|
|
|
|
!
|
2005-12-19 21:31:56 +08:00
|
|
|
USE parameters, ONLY : natx
|
|
|
|
USE kinds, ONLY : DP
|
|
|
|
!
|
|
|
|
#ifdef DFT_CP
|
2005-11-25 05:17:29 +08:00
|
|
|
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
|
2005-12-14 20:05:37 +08:00
|
|
|
USE wave_base, ONLY : frice
|
2005-11-25 05:17:29 +08:00
|
|
|
USE energies, ONLY : eself, etot
|
|
|
|
USE from_scratch_module, ONLY : from_scratch
|
|
|
|
USE gvecs, ONLY : ngs
|
|
|
|
USE ions_positions, ONLY : tau0
|
2005-12-19 21:31:56 +08:00
|
|
|
USE ions_base, ONLY : ityp, nat
|
2005-11-25 05:17:29 +08:00
|
|
|
USE phase_factors_module, ONLY : strucf
|
|
|
|
USE reciprocal_vectors, ONLY : mill_l
|
|
|
|
USE time_step, ONLY : dt2
|
2005-12-19 21:31:56 +08:00
|
|
|
USE io_global, ONLY : stdout
|
|
|
|
#endif
|
|
|
|
!
|
|
|
|
#ifdef DFT_PW
|
|
|
|
USE force_mod, ONLY : force
|
|
|
|
USE ener, ONLY : etot
|
|
|
|
USE cell_base, ONLY : alat
|
|
|
|
USE constants, ONLY : e2
|
|
|
|
#endif
|
2005-11-25 05:17:29 +08:00
|
|
|
IMPLICIT NONE
|
|
|
|
!
|
|
|
|
! ... input variables
|
|
|
|
!
|
|
|
|
REAL (KIND=DP) :: fion(3,natx)
|
|
|
|
!
|
2005-12-19 21:31:56 +08:00
|
|
|
#ifdef DFT_CP
|
|
|
|
!
|
2005-11-25 05:17:29 +08:00
|
|
|
! ... local variables
|
|
|
|
!
|
|
|
|
LOGICAL :: tfirst, tlast
|
2005-12-14 20:05:37 +08:00
|
|
|
REAL (KIND=DP) :: enthal, fccc
|
2005-11-25 05:17:29 +08:00
|
|
|
REAL (KIND=DP) :: dt2bye, enb, enbi, ccc
|
2005-12-19 21:31:56 +08:00
|
|
|
integer :: ipol, na
|
2005-11-25 05:17:29 +08:00
|
|
|
!
|
|
|
|
! ... for smooth restart in the new coordinates
|
|
|
|
!
|
2005-12-19 21:31:56 +08:00
|
|
|
fion = 0.0
|
2005-11-25 05:17:29 +08:00
|
|
|
dt2bye = dt2 / emass
|
2005-12-14 20:05:37 +08:00
|
|
|
fccc = 1.D0 / ( 1.D0 + frice )
|
2005-11-25 05:17:29 +08:00
|
|
|
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 )
|
|
|
|
!
|
2005-12-14 20:05:37 +08:00
|
|
|
! ... relax wavefunction in new position
|
2005-11-25 05:17:29 +08:00
|
|
|
!
|
|
|
|
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
|
2005-12-19 21:31:56 +08:00
|
|
|
!
|
|
|
|
! ... write on output the forces
|
|
|
|
!
|
|
|
|
DO na = 1, nat
|
|
|
|
WRITE( stdout, 9035) na, ityp(na), ( fion(ipol,na), ipol = 1, 3 )
|
|
|
|
END DO
|
|
|
|
9035 FORMAT(5X,'atom ',I3,' type ',I2,' force = ',3F14.8)
|
|
|
|
|
|
|
|
!
|
|
|
|
#endif
|
|
|
|
!
|
|
|
|
#ifdef DFT_PW
|
|
|
|
fion = 0.0
|
|
|
|
force = 0.0
|
2005-12-25 07:14:27 +08:00
|
|
|
!call hinit0 ()
|
2005-12-19 21:31:56 +08:00
|
|
|
CALL hinit1 ()
|
|
|
|
call electrons()
|
|
|
|
call forces()
|
|
|
|
! convert to hartree/bohr
|
|
|
|
fion = force /e2
|
|
|
|
etot = etot / e2
|
|
|
|
#endif
|
2005-11-25 05:17:29 +08:00
|
|
|
!
|
|
|
|
RETURN
|
|
|
|
END SUBROUTINE relax_wavefunction
|
2005-12-23 07:35:20 +08:00
|
|
|
!
|
|
|
|
! ----------------------------------------------
|
|
|
|
!
|
2005-12-19 21:31:56 +08:00
|
|
|
subroutine set_guess_wfc ( disp_sign )
|
|
|
|
!
|
|
|
|
USE kinds, ONLY : DP
|
|
|
|
#ifdef DFT_CP
|
|
|
|
USE ions_positions, ONLY : tau0
|
|
|
|
USE cp_main_variables, ONLY : lambda, lambdam, nfi
|
|
|
|
USE wavefunctions_module, ONLY : c0, cm
|
|
|
|
USE vibrations, ONLY : ref_c0
|
|
|
|
#endif
|
|
|
|
#ifdef DFT_PW
|
|
|
|
USE scf, ONLY : rho
|
|
|
|
USE wavefunctions_module, ONLY : evc
|
|
|
|
USE io_files, ONLY : nwordwfc, iunwfc, iunoldwfc2, prefix
|
2005-12-25 07:14:27 +08:00
|
|
|
USE ions_base, ONLY : nat, ntyp => nsp, ityp, tau
|
2005-12-19 21:31:56 +08:00
|
|
|
USE klist, ONLY : nks
|
2005-12-25 07:14:27 +08:00
|
|
|
USE scf, ONLY : rho, rho_core, vr
|
|
|
|
USE gvect, ONLY : nrxx, ngm, g, gg, gstart, &
|
|
|
|
nr1, nr2, nr3, nl, &
|
|
|
|
eigts1, eigts2, eigts3, &
|
|
|
|
nrx1, nrx2, nrx3
|
|
|
|
USE cell_base, ONLY : omega, bg, alat
|
|
|
|
USE ener, ONLY : ehart, etxc, vtxc
|
|
|
|
USE extfield, ONLY : etotefield
|
|
|
|
USE vlocal, ONLY : strf
|
2005-12-19 21:31:56 +08:00
|
|
|
|
|
|
|
#endif
|
|
|
|
!
|
|
|
|
! ... input variables
|
|
|
|
!
|
|
|
|
INTEGER, INTENT(IN) :: disp_sign
|
|
|
|
!
|
|
|
|
! ... local variables
|
|
|
|
!
|
|
|
|
logical :: exst
|
2005-12-25 07:14:27 +08:00
|
|
|
REAL(kind=DP) :: charge
|
2005-12-19 21:31:56 +08:00
|
|
|
!
|
|
|
|
!
|
|
|
|
#ifdef DFT_CP
|
|
|
|
nfi = 0
|
|
|
|
#endif
|
|
|
|
!
|
|
|
|
!
|
|
|
|
|
|
|
|
IF(disp_sign.EQ.-1) THEN
|
|
|
|
!
|
|
|
|
! ... use reference wavefunction as initial guess
|
|
|
|
!
|
|
|
|
#ifdef DFT_CP
|
|
|
|
c0(:,:,1,1) = ref_c0(:,:,1,1)
|
|
|
|
#endif
|
|
|
|
!
|
|
|
|
#ifdef DFT_PW
|
|
|
|
INQUIRE (file=TRIM(prefix)//'.old2rho', EXIST=exst)
|
|
|
|
if ( exst ) then
|
|
|
|
!
|
|
|
|
! ... charge-density from file to memory
|
|
|
|
!
|
|
|
|
CALL io_pot( - 1, 'old2rho', rho, 1 )
|
|
|
|
!
|
|
|
|
CALL v_of_rho( rho, rho_core, nr1, nr2, nr3, nrx1, nrx2, nrx3, &
|
2005-12-23 07:35:20 +08:00
|
|
|
nrxx, nl, ngm, gstart, nspin, g, gg, alat, omega, &
|
|
|
|
ehart, etxc, vtxc, etotefield, charge, vr )
|
2005-12-19 21:31:56 +08:00
|
|
|
!
|
|
|
|
! .. swapping the wavefunctions
|
|
|
|
!
|
|
|
|
CALL diropn( iunoldwfc2, 'old2wfc', nwordwfc, exst )
|
|
|
|
!
|
|
|
|
DO ik = 1, nks
|
|
|
|
!
|
|
|
|
! ... "2old" -> "now"
|
|
|
|
!
|
|
|
|
CALL davcio( evc, nwordwfc, iunoldwfc2, ik, - 1 )
|
|
|
|
CALL davcio( evc, nwordwfc, iunwfc, ik, + 1 )
|
|
|
|
!
|
|
|
|
END DO
|
|
|
|
!
|
|
|
|
CLOSE( UNIT = iunoldwfc2 )
|
|
|
|
!
|
|
|
|
end if
|
|
|
|
history = 1
|
|
|
|
call update_pot()
|
|
|
|
!
|
|
|
|
#endif
|
|
|
|
!
|
|
|
|
ELSE
|
|
|
|
!
|
|
|
|
! ... extrapolate wave function coefficients from their value
|
|
|
|
! ... at dispalcement at the negative direction
|
|
|
|
!
|
|
|
|
#ifdef DFT_CP
|
|
|
|
c0(:,:,1,1) = 2*ref_c0(:,:,1,1)-c0(:,:,1,1)
|
|
|
|
#endif
|
|
|
|
#ifdef DFT_PW
|
|
|
|
history = 2
|
|
|
|
call update_pot()
|
|
|
|
#endif
|
|
|
|
!
|
|
|
|
END IF
|
|
|
|
!
|
|
|
|
! ... set wavefunction velocity to zero
|
|
|
|
!
|
|
|
|
#ifdef DFT_CP
|
|
|
|
cm = c0
|
|
|
|
lambdam = lambda
|
|
|
|
#endif
|
2006-02-08 20:42:32 +08:00
|
|
|
!
|
2005-12-19 21:31:56 +08:00
|
|
|
return
|
|
|
|
end subroutine set_guess_wfc
|
2005-12-23 07:35:20 +08:00
|
|
|
!
|
|
|
|
!
|
|
|
|
! ----------------------------------------------
|
|
|
|
!
|
|
|
|
!
|
|
|
|
SUBROUTINE cofmass( tau, mass, nat, com )
|
|
|
|
!
|
|
|
|
! ... calculation of the center of mass
|
|
|
|
!
|
|
|
|
USE kinds, only : DP
|
|
|
|
IMPLICIT NONE
|
|
|
|
REAL(DP), INTENT(IN) :: tau(3,nat), mass(nat)
|
|
|
|
REAL(DP), INTENT(OUT) :: com(3)
|
|
|
|
INTEGER, INTENT(IN) :: nat
|
|
|
|
|
|
|
|
REAL(DP) :: tmas
|
|
|
|
INTEGER :: is, i, ia
|
|
|
|
!
|
|
|
|
tmas=SUM(mass(1:nat))
|
|
|
|
!
|
|
|
|
do i=1,3
|
|
|
|
com(i)=DOT_PRODUCT(tau(i,1:nat),mass(1:nat))
|
|
|
|
com(i)=com(i)/tmas
|
|
|
|
end do
|
|
|
|
!
|
|
|
|
RETURN
|
|
|
|
END SUBROUTINE cofmass
|
2005-12-29 16:38:27 +08:00
|
|
|
|
|
|
|
!-----------------------------------------------------------------------
|
|
|
|
SUBROUTINE vib_pbc(rin,a1,a2,a3,ainv,rout)
|
|
|
|
!-----------------------------------------------------------------------
|
|
|
|
!
|
|
|
|
! brings atoms inside the unit cell
|
|
|
|
!
|
|
|
|
IMPLICIT NONE
|
|
|
|
! input
|
|
|
|
REAL(8) rin(3), a1(3),a2(3),a3(3), ainv(3,3)
|
|
|
|
! output
|
|
|
|
REAL(8) rout(3)
|
|
|
|
! local
|
|
|
|
REAL(8) x,y,z
|
|
|
|
!
|
|
|
|
! bring atomic positions to crystal axis
|
|
|
|
!
|
|
|
|
x = ainv(1,1)*rin(1)+ainv(1,2)*rin(2)+ainv(1,3)*rin(3)
|
|
|
|
y = ainv(2,1)*rin(1)+ainv(2,2)*rin(2)+ainv(2,3)*rin(3)
|
|
|
|
z = ainv(3,1)*rin(1)+ainv(3,2)*rin(2)+ainv(3,3)*rin(3)
|
|
|
|
!
|
|
|
|
! bring x,y,z in the range between -0.5 and 0.5
|
|
|
|
!
|
|
|
|
x = x - NINT(x)
|
|
|
|
y = y - NINT(y)
|
|
|
|
z = z - NINT(z)
|
|
|
|
!
|
|
|
|
! bring atomic positions back in cartesian axis
|
|
|
|
!
|
|
|
|
rout(1) = x*a1(1)+y*a2(1)+z*a3(1)
|
|
|
|
rout(2) = x*a1(2)+y*a2(2)+z*a3(2)
|
|
|
|
rout(3) = x*a1(3)+y*a2(3)+z*a3(3)
|
|
|
|
!
|
|
|
|
RETURN
|
|
|
|
END SUBROUTINE vib_pbc
|