!======================================================================= ! program cpe ! !======================================================================= !*** Molecular Dynamics using Density-Functional Theory **** !*** this is a Car-Parrinello program using Vanderbilt pseudopotentials !*********************************************************************** !*** based on version 11 of cpv code including ggapw 07/04/99 !*** copyright Alfredo Pasquarello 10/04/1996 !*** parallelized and converted to f90 by Paolo Giannozzi (2000), !*** using parallel FFT written for PWSCF by Stefano de Gironcoli !*** PBE added by Michele Lazzeri (2000) !*** finite electric field Paolo Umari(2001) !*********************************************************************** !*** appropriate citation for use of this code: !*** Car-Parrinello method R. Car and M. Parrinello PRL 55, 2471 (1985) !*** current implementation A. Pasquarello, K. Laasonen, R. Car, !*** C. Lee, and D. Vanderbilt, PRL 69, 1982 (1992); !*** K. Laasonen, A. Pasquarello, R. Car, !*** C. Lee, and D. Vanderbilt, PRB 47, 10142 (1993). !*** implementation gga A. Dal Corso, A. Pasquarello, A. Baldereschi, !*** and R. Car, PRB 53, 1180 (1996). !*********************************************************************** !*** !*** f90 version, with dynamical allocation of memory !*** Variables that do not change during the dynamics are in modules !*** All other variables are passed as arguments !*** !*********************************************************************** !*** !*** fft : uses machine's own complex fft routines, two real fft at the time !*** ggen: g's only for positive halfspace (g>) !*** all routines : keep equal c(g) and [c(-g)]* !*** !*********************************************************************** ! general variables: ! delt = delta t ! emass = electron mass (fictitious) ! dt2bye = 2*delt/emass !*********************************************************************** ! ! ATTENZIONE l'iforce diventa vettore per tutte le componenti cartesiane use control use core use cnst use cvan use ener use elct use gvec use gvecs use gvecb use ions use parm use parms use parmb use pseu use work1 use work_box use work2 #ifdef PARA use para_mod use work_fft #endif implicit none ! ! control variables ! logical trane,twall,tbump,tnosee,tortho logical tfor,tnosep,trhor,trhow logical tfirst,tbuff,tlast,tvlocw ! ! wavefunctions ! complex(kind=8), allocatable:: c0(:,:), cm(:,:), phi(:,:) ! ! structure factors e^{-ig*R} ! complex(kind=8), allocatable:: ei1(:,:,:), ei2(:,:,:), ei3(:,:,:) complex(kind=8), allocatable:: eigr(:,:,:) ! ! structure factors (summed over atoms of the same kind) ! complex(kind=8), allocatable:: sfac(:,:) ! ! indexes, positions, and structure factors for the box grid ! integer irb(3,nax,nsx) real(kind=8) taub(3,nax,nsx) complex(kind=8), allocatable:: eigrb(:,:,:) ! ! charge densities and potentials ! rhog = charge density in g space ! rhor = charge density in r space (dense grid) ! rhos = charge density in r space (smooth grid) ! complex(kind=8), allocatable:: rhog(:,:) real(kind=8), allocatable:: rhor(:,:), rhos(:,:) ! ! nonlocal projectors: ! bec = scalar product of projectors and wave functions ! betae = nonlocal projectors in g space = beta x e^(-ig.R) ! becdr = used in force calculation ! complex(kind=8), allocatable:: betae(:,:) real(kind=8), allocatable:: bec(:,:), becdr(:,:,:) real(kind=8), allocatable:: bephi(:,:), becp(:,:) real(kind=8), allocatable:: rhovan(:,:,:,:), deeq(:,:,:,:) ! ! mass preconditioning ! real(kind=8), allocatable:: ema0bg(:) real(kind=8) emaec ! ! constraints (lambda at t, lambdam at t-dt, lambdap at t+dt) ! real(kind=8), allocatable:: lambda(:,:), lambdam(:,:), lambdap(:,:) ! ! ionic positions, center of mass position ! real(kind=8) tau0(3,nax,nsx), taum(3,nax,nsx), taup(3,nax,nsx), & & taui(3,nax,nsx) real(kind=8) cdmp(3),cdm0(3),cdmm(3),cdmvel(3) ! ! forces on ions ! real(kind=8) fion(3,nax,nsx), fionm(3,nax,nsx) integer iforce(3,nax,nsx) ! real(kind=8) vel(3,nax,nsx) ! integer ndr, ndw, nbeg, maxit, nomore, ntop, iprint ! ! variables for Berry phase ! integer, allocatable:: ctable(:,:,:)!correspondence tables for diff. polarization integer, allocatable:: ctabin(:,:,:)!inverse correspondence table complex(kind=8), allocatable:: qmat(:,:)!inverse of matrix Q, for Barry's phase complex(kind=8), allocatable:: gqq(:,:,:,:)!factors int beta_Ri^*beta_Rj exp(iGr)dr complex(kind=8), allocatable:: gqqm(:,:,:,:)! the same with exp(-iGr) complex(kind=8), allocatable:: gqq0(:,:,:,:)!factors int beta_Ri^*beta_Rj exp(iGr)dr, at Gamma complex(kind=8), allocatable:: gqqm0(:,:,:,:)! the same with exp(-iGr), at Gamma real(kind=8) evalue!strenght of electric field integer ipol!direction of electric field complex(kind=8), allocatable:: df(:) complex(kind=8) detq real(kind=8) enb!berry energy, electronic part real(kind=8) enbi!berry energy, ionic part integer unit! unit for stdin integer istat!error value for flush real(kind=8) cdzp(3),cdzm(3), cdz0(3)!centers of ionic charges real(kind=8) taunew(3,nax,nsx)!the electric field must be always along the column direction !this just in the parallel case integer ipolp!this funny parameter is used to impose z-direction in parallel case real(kind=8) pberryel, pberryion!polarizations ! ! work variables ! real(kind=8) acc(10) complex(kind=8), allocatable:: c2(:), c3(:) complex(kind=8) speed real(kind=8) & & tempp, xnhe0, vnhp, xnhp0, xnhpm, verl1, verl2, verl3, & & velrel, fccc, xnhem, vnhe, anor, savee, savep, press, & & enthal, epot, xnhpp, xnhep, epre, enow, tps, econs, econt, & & ampre, fricp, greasp, eps, reps, qnp, tempw, qne, factem, & & frice, grease, emass, delt, ccc, bigr, fordt2, dt2, & & dt2by2, twodel, gausp, dt2bye, volm, vol0, gkbt, dt2hbe, & & dtbyq, dt2byq real(kind=8) ekinv, ekinc0, ekinp, ekinpr, ekincm, ekinc, ekincw integer nnn, is, nacc, ia, j, iter, nfi, ig, nr, i, iv,jv, & & isa, nerre, iss ,flag ! ! program starts here ! #if defined(T3E) || defined(ORIGIN) || defined(OSF1) unit=9 #else unit=5 #endif #ifdef PARA call startup #endif ! ! ================================================================== ! ==== 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 ==== ! ================================================================== scmass = 1822.89d0 factem = 27.212d0*11605.d0 ! ! ================================================================== ! read control variables from stdin (unit 5) ! ================================================================== call init0 (nbeg,ndr,ndw,nomore,iprint & & ,delt,emass,emaec & & ,frice,grease,twall & & ,tortho,eps,maxit,reps,ntop & & ,trane,ampre & & ,tfor,fricp,greasp & & ,trhor,trhow,tvlocw & & ,tnosep,qnp,tempw & & ,tnosee,qne,ekincw) ! ================================================================== ! ! general variables ! tfirst=.true. nacc=5 ! gausp=delt*sqrt(tempw/factem) twodel=2.d0*delt fordt2=twodel*twodel dt2=delt*delt dt2by2=.5d0*dt2 dt2bye=dt2/emass dt2hbe=dt2by2/emass dtbyq=delt/qnp dt2byq=dt2/qnp #ifdef PARA call openfiles(tbuff,nbeg,ndr,ndw) #endif ! write(6,*) ' out from init0' ! ================================================================== ! read structural data (cell, atomic positions, etc) ! initialize g-vectors, fft grids ! ================================================================== call init1 (taunew,tbuff,iforce) read(unit,*) ipol, evalue write(6,*) 'ipol', ipol, 'Electric Field', evalue ! ================================================================== ! rearrange carthesian axes ! ================================================================== #ifdef PARA do is=1,nsp do ia=1,na(is) if(ipol .eq. 1) then tau0(1,ia,is)=taunew(2,ia,is) tau0(2,ia,is)=taunew(3,ia,is) tau0(3,ia,is)=taunew(1,ia,is) elseif( ipol .eq. 2) then tau0(1,ia,is)=taunew(3,ia,is) tau0(2,ia,is)=taunew(1,ia,is) tau0(3,ia,is)=taunew(2,ia,is) elseif( ipol .eq. 3) then tau0(1,ia,is)=taunew(1,ia,is) tau0(2,ia,is)=taunew(2,ia,is) tau0(3,ia,is)=taunew(3,ia,is) end if end do end do ipolp=3 #else do is=1,nsp do ia=1,na(is) tau0(1,ia,is)=taunew(1,ia,is) tau0(2,ia,is)=taunew(2,ia,is) tau0(3,ia,is)=taunew(3,ia,is) end do end do ipolp=ipol #endif ! write(6,*) ' out from init1' !----------read ipol and evalue from stdin------------- ! ! more initialization requiring atomic positions ! nat=0 nas=0 do is=1,nsp nas=max(nas,na(is)) nat=nat+na(is) end do call cofmass(tau0,cdm0) if(iprsta.gt.1) then write(6,*) ' tau0 ' write(6,*) (((tau0(i,ia,is),i=1,3),ia=1,na(is)),is=1,nsp) endif do i=1,3 cdmm(i)=cdm0(i) cdmvel(i)=0. end do do is=1,nsp do ia=1,na(is) do i=1,3 taum(i,ia,is)=tau0(i,ia,is) taup(i,ia,is)=0.0 end do end do end do volm=omega vol0=omega gkbt = 3.*nat*tempw/factem ! ================================================================== ! allocate and initialize nonlocal potentials ! ================================================================== call nlinit write(6,*) ' out from nlinit' ! ! ================================================================== ! allocation of all arrays not already allocated in init1 and nlinit ! ================================================================== ! allocate(c0(ngw,nx)) allocate(cm(ngw,nx)) allocate(phi(ngw,nx)) write(6,*) 'allo1' allocate(wrk2(ngw,max(nas,n))) allocate(eigr(ngw,nas,nsp)) allocate(eigrb(ngb,nas,nsp)) allocate(sfac(ngs,nsp)) allocate(rhops(ngs,nsp)) allocate(vps(ngs,nsp)) allocate(rhor(nnr,nspin)) allocate(rhos(nnrs,nspin)) allocate(rhog(ng,nspin)) write(6,*) 'allo2' allocate(wrk1(nnr)) allocate(qv(nnrb)) allocate(c2(ngw)) allocate(c3(ngw)) allocate(ema0bg(ngw)) allocate(lambda(nx,nx)) allocate(lambdam(nx,nx)) allocate(lambdap(nx,nx)) write(6,*) 'allo3' allocate(ei1(-nr1:nr1,nas,nsp)) allocate(ei2(-nr2:nr2,nas,nsp)) allocate(ei3(-nr3:nr3,nas,nsp)) write(6,*) 'allo4' write(6,*) ngw,nhsa,n,nat,nhx,nspin allocate(betae(ngw,nhsa)) allocate(becdr(nhsa,n,3)) allocate(bec (nhsa,n)) allocate(bephi(nhsa,n)) allocate(becp (nhsa,n)) allocate(deeq(nat,nhx,nhx,nspin)) allocate(rhovan(nat,nhx,nhx,nspin)) ! berry phase arrays write(6,*) 'allo5' write(6,*) ngw,nx,nhx,nas,nsp allocate( ctable(ngw,2,3)) allocate( ctabin(ngw,2,3)) allocate( qmat(nx,nx)) allocate( gqq(nhx,nhx,nas,nsp)) allocate( gqqm(nhx,nhx,nas,nsp)) allocate( df(ngw)) allocate( gqq0(nhx,nhx,nas,nsp)) allocate( gqqm0(nhx,nhx,nas,nsp)) write(6,*) 'allo6' #ifdef PARA allocate(aux(nnr)) #endif call zero(nat*nhx*nhx*nspin,deeq) do i=1,n call zero(n,lambda(1,i)) call zero(2*ngw,cm(1,i)) call zero(2*ngw,c0(1,i)) end do ! ! local potential initialization ! call formf ! ! mass preconditioning: ema0bg(i) = ratio of emass(g=0) to emass(g) ! for g**2>emaec the electron mass ema0bg(g) rises quadratically ! do i=1,ngw ema0bg(i)=1./max(1.d0,tpiba2*g(i)/emaec) if(iprsta.ge.10) write(6,*) i,' ema0bg(i) ',ema0bg(i) end do ! if (nbeg.lt.0) then !======================================================================= ! nbeg = -1 or nbeg = -2 or nbeg = -3 !======================================================================= ! !======================================================================= ! nbeg = -1 or nbeg = -2 or nbeg = -3 !======================================================================= ! if(nbeg.eq.-1) then flag=0 ! this is the value for normal operation ! flag=1 ! only for LDOS !!! AAAAAA #ifdef PARA call readpfile & #else call readfile & #endif & ( flag,nx,nax,nsp,ndr,nfi,ngw,n,nr,cm,cm,tau0,vol0,cdm0,cdmm, & & taum,volm,acc,taui,a1,a2,a3,alat,lambda,lambdam, & & xnhe0,xnhem,vnhe,xnhp0,xnhpm,vnhp,vel,ekincm) endif ! ! phfac calculates the factors ei1, ei2, ei3, and eigr=ei1*ei2*ei3 ! call phfac(tau0,ei1,ei2,ei3,eigr) ! call initbox ( nr1, nr2, nr3, nr1b, nr2b, nr3b, nsp, na, & & a1, a2, a3, ainv, tau0, taub, irb ) ! call phbox(taub,eigrb) ! !################################################################## ! the flag above must be modified accordingly !! ! ! call eigrot(nspin,nx,nupdwn,iupdwn,f,lambda,c0,cm) ! call calbec (1,nsp,eigr,c0,bec) ! call eigfun (nfi,iprint,tbuff,trhor,trhow,& ! &c0,irb,eigrb,bec,rhovan,rhor,rhog,rhos,& ! &2,1,350) !#ifdef PARA ! call mpi_finalize(i) !#endif ! stop !################################################################## ! if(iprsta.gt.2)then do is=1,nvb write(6,'(/,2x,''species= '',i2)') is do ia=1,na(is) write(6,2000) ia, (irb(i,ia,is),i=1,3) 2000 format(2x,'atom= ',i3,' irb1= ',i3,' irb2= ',i3, & & ' irb3= ',i3) end do end do endif ! if(trane) then ! ! random initialization ! call COPY(2*ngw*n,cm,1,c0,1) call randin(1,n,ng0,ngw,ampre,cm) else if (nbeg.eq.-3) then ! ! gaussian initialization ! call gausin(eigr,cm) endif ! ! prefor calculates betae (used by graham) ! call prefor(eigr,betae) call graham(betae,bec,cm) if(iprsta.ge.3) call dotcsc(eigr,cm) ! nfi=0 ! ! strucf calculates the structure factor sfac ! call strucf (ei1,ei2,ei3,sfac) call calbec (nvb+1,nsp,eigr,cm,bec) ! call rhoofr (nfi,iprint,tbuff,trhor,trhow, & & cm,irb,eigrb,bec,rhovan,rhor,rhog,rhos) if(iprsta.gt.0) write(6,*) ' out from rhoofr' call vofrho(nfi,iprint,rhor,rhog,rhos,tfor,tfirst,tlast, & & ei1,ei2,ei3,sfac,tau0,tvlocw,fion) if(iprsta.gt.0) write(6,*) ' out from vofrho' if(iprsta.gt.2) write(6,*) ' fion ', & & (((fion(i,ia,is),i=1,3),ia=1,na(is)),is=1,nsp) ! set up of Berry stuff write(6,'(''before gtable'')') call gtable(ipolp,ctable(1,1,ipolp)) write(6,'(''out of gtable'')') call gtablein(ipolp,ctabin(1,1,ipolp)) write(6,'(''out of gtablein'')') call qqberry2(gqq0,gqqm0,ipolp)!for Vanderbilt pps write(6,'(''out of qqberry2'')') call qqupdate(eigr,gqqm0,gqq,gqqm,ipolp) write(6,'(''out of qqupdate'')') ! computes Q^-1 and detQ call qmatrixd(cm,bec,ctable(1,1,ipolp),gqq,qmat,detq) write(6,'(''out of qmatrixd'')') call enberry( detq, ipolp,enb) call berryion(tau0,fion,tfor,ipolp,evalue,enbi) pberryel=enb pberryion=enbi enb=enb*evalue enbi=enbi*evalue etot=etot+enb+enbi if(mod(nfi-1,iprint).eq.0 .or. tlast) & & write(6,1956) 'Polarizzazione :', nfi, pberryel,pberryion,pberryel+pberryion ! ! forces for eigenfunctions ! ! newd calculates deeq and a contribution to fion ! call newd(rhor,tfor,irb,eigrb,rhovan,deeq,fion) write(6,*) ' out from newd' call prefor(eigr,betae) ! ! if n is odd => c(*,n+1)=0 ! do i=1,n,2 call dforce & & (bec,deeq,betae,i,cm(1,i),cm(1,i+1),c2,c3,rhos,tbuff) call dforceb & & (cm, i, betae, ipolp, bec ,ctabin(1,1,ipolp), gqq, gqqm, qmat, deeq, df) do ig=1,ngw c2(ig)=c2(ig)+evalue*df(ig) enddo call dforceb & & (cm, i+1, betae, ipolp, bec ,ctabin(1,1,ipolp), gqq, gqqm, qmat, deeq, df) do ig=1,ngw c3(ig)=c3(ig)+evalue*df(ig) enddo ccc=dt2hbe do j=1,ngw c0(j,i) =cm(j,i) +ccc*ema0bg(j)*c2(j) c0(j,i+1)=cm(j,i+1)+ccc*ema0bg(j)*c3(j) end do end do write(6,*) ' out from dforce' ! ! buffer for wavefunctions is unit 21 ! if(tbuff) rewind 21 ! ! nlfq needs deeq calculated in newd ! if (tfor) call nlfq(cm,deeq,eigr,bec,becdr,fion) write(6,*) ' out from nlfq' ! ========================================================== ! imposing the orthogonality ! ========================================================== ! call calphi(cm,ema0bg,bec,betae,phi) write(6,*) ' out from calphi' ! if (tortho) then call ortho (eigr,c0,phi,lambda, & & bigr,iter,ccc,eps,maxit,delt,bephi,becp) else call graham(betae,bec,c0) write(6,*) ' graham c0 ' end if ! ! nlfl needs lambda becdr and bec ! if (tfor) call nlfl(bec,becdr,lambda,fion) write(6,*) ' out from nlfl' ! ! bforceion adds the force term due to electronic berry phase ! only in US-case call bforceion(fion,tfor,ipolp, qmat,bec,becdr,gqq,evalue) write(6,*) ' out from bforceion' if(iprsta.ge.3) then nnn=min(12,n) write(6,*) 'from main:' do i=1,nnn write(6,'(12f8.5)') (lambda(i,j)*ccc,j=1,nnn) end do write(6,*) endif ! if (tortho) then call updatc(ccc,lambda,phi,bephi,becp,bec,c0) write(6,*) ' out from updatc' end if call calbec (nvb+1,nsp,eigr,c0,bec) !------------------------------------------------ ! the center of mass along ipolp is kept steady !------------------------------------------------ !if(tfor) call noforce(fion, ipolp) ! write(6,*) ' out from calbec' ! ================================================================== ! cm now orthogonalized ! ================================================================== if(iprsta.ge.3) call dotcsc(eigr,c0) ! if(tfor) then do is=1,nsp do ia=1,na(is) do i=1,3 tau0(i,ia,is)=taum(i,ia,is) & & +iforce(i,ia,is)*dt2by2 & & /pmass(is)*fion(i,ia,is) end do end do end do ! ==================================================================== ! the center of ionic charge is kept steady along direction ipolp ! ==================================================================== ! call cofcharge(tau0,cdzp) ! call cofcharge(taum,cdz0) ! do is=1,nsp ! do ia=1,na(is) ! do i=1,3 ! tau0(i,ia,is)=tau0(i,ia,is)-cdzp(i)+cdz0(i) ! enddo ! enddo ! enddo ! call keptwfc(c0,-cdzp(1)+cdz0(1),-cdzp(2)+cdz0(2),-cdzp(3)+cdz0(3)) ! call phfac(tau0,ei1,ei2,ei3,eigr) ! !call calbec (1,nsp,eigr,c0,bec) endif ! if (tnosep) then xnhp0=0. xnhpm=0. vnhp =0. do is=1,nsp do ia=1,na(is) do i=1,3 fionm(i,ia,is)=0. vel(i,ia,is)=(tau0(i,ia,is)-taum(i,ia,is))/delt end do end do end do endif ! ! ====================================================== ! kinetic energy of the electrons ! ====================================================== ekincm=0.0 do i=1,n do j=1,ngw speed=c0(j,i)-cm(j,i) ekincm=ekincm+2.*real(conjg(speed)*speed)/ema0bg(j) end do if (ng0.eq.2) then speed=c0(1,i)-cm(1,i) ekincm=ekincm-real(conjg(speed)*speed) end if end do ekincm=ekincm*emass/dt2 #ifdef PARA call reduce(1,ekincm) #endif ! if(tnosee)then xnhe0=0. xnhem=0. vnhe =0. endif ! do i=1,n do j=1,n lambdam(j,i)=lambda(j,i) end do end do ! else !======================================================================= ! nbeg = 0 or nbeg = 1 !======================================================================= #ifdef PARA call readpfile & #else call readfile & #endif & ( 1,nx,nax,nsp,ndr,nfi,ngw,n,nr,c0,cm,tau0,vol0,cdm0,cdmm, & & taum,volm,acc,taui,a1,a2,a3,alat,lambda,lambdam, & & xnhe0,xnhem,vnhe,xnhp0,xnhpm,vnhp,vel,ekincm) write(6,*) ' disk read ',ndr ! call phfac(tau0,ei1,ei2,ei3,eigr) if(trane.and.trhor) then call randin(nr+1,n,ng0,ngw,ampre,c0) call prefor(eigr,betae) call graham(betae,bec,c0) call COPY(2*ngw*n,c0,1,cm,1) endif ! if(iprsta.gt.2) then write(6,*) ' read: tau0 ' write(6,*) (((tau0(i,ia,is),i=1,3),ia=1,na(is)),is=1,nsp) endif ! omega=vol0 tpiba=2.*pi/alat tpiba2=tpiba*tpiba ! call strucf (ei1,ei2,ei3,sfac) call calbec (1,nsp,eigr,c0,bec) ! set up of Berry stuff call gtable(ipolp,ctable(1,1,ipolp)) call gtablein(ipolp,ctabin(1,1,ipolp)) write(6,'(''out of gtable'')') call qqberry2(gqq0,gqqm0,ipolp)!for Vanderbilt pps call qqupdate(eigr,gqqm0,gqq,gqqm,ipolp) call cofcharge(tau0,cdz0) ! end if !==============================================end of if(nbeg.lt.0)===== ! ! ================================================================== ! restart with new averages and nfi=0 ! ================================================================== if(nbeg.le.0) then call cofmass(tau0,cdm0) do i=1,10 acc(i)=0. end do do is=1,nsp do ia=1,na(is) do i=1,3 taui(i,ia,is)=tau0(i,ia,is) end do end do end do nfi=0 end if ! if(.not.tfor) then do is=1,nsp do ia=1,na(is) do i=1,3 fion(i,ia,is)=0.d0 end do end do end do end if ! fccc=1. nomore=nomore+nfi ! !====================================================================== ! ! basic loop for molecular dynamics starts here ! !====================================================================== ! 1000 continue call tictac(1,0) ! ! calculation of velocity of nose-hoover variables ! if(tnosep)then vnhp=2.*(xnhp0-xnhpm)/delt-vnhp do is=1,nsp do ia=1,na(is) do i=1,3 vel(i,ia,is) = 2.*(tau0(i,ia,is)-taum(i,ia,is))/delt & & - vel(i,ia,is) end do end do end do endif fccc=1./(1.+frice) if(tnosee)then vnhe=2.*(xnhe0-xnhem)/delt-vnhe fccc=1./(1.+0.5*delt*vnhe) endif ! if (tfor.or.tfirst) then call initbox ( nr1, nr2, nr3, nr1b, nr2b, nr3b, nsp, na, & & a1, a2, a3, ainv, tau0, taub, irb ) call phbox(taub,eigrb) endif ! if(tfor) call phfac(tau0,ei1,ei2,ei3,eigr) ! ! strucf calculates the structure factor sfac ! call strucf (ei1,ei2,ei3,sfac) ! nfi=nfi+1 tlast=(nfi.eq.nomore) ! call rhoofr (nfi,iprint,tbuff,trhor,trhow, & & c0,irb,eigrb,bec,rhovan,rhor,rhog,rhos) ! ! charge density, must appear right after the rhoofr call ! call totrho(rhor) ! electric field along c-axis call efield(rhog) ! call vofrho(nfi,iprint,rhor,rhog,rhos,tfor,tfirst,tlast, & & ei1,ei2,ei3,sfac,tau0,tvlocw,fion) ! average total potential in r-space ! call totpot(rhor) if( tfor ) call qqupdate(eigr,gqqm0,gqq,gqqm,ipolp) ! computes Q^-1 and detQ if((evalue .ne. 0.d0) .or. (mod(nfi-1,10).eq.0.) .or. tfirst) then call qmatrixd(c0,bec,ctable(1,1,ipolp),gqq,qmat,detq) call enberry( detq, ipolp,enb) call berryion(tau0,fion,tfor,ipolp,evalue,enbi) pberryel=enb pberryion=enbi endif enb=enb*evalue enbi=enbi*evalue etot=etot+enb+enbi if(mod(nfi-1,iprint).eq.0 .or. tlast) & & write(6,1956) 'Polarizzazione :', nfi, pberryel,pberryion,pberryel+pberryion write(70,*) nfi,pberryel+pberryion !###################################################################### ! call eigrot(nspin,nx,nupdwn,iupdwn,f,lambda,c0,cm) ! call calbec (1,nsp,eigr,c0,bec) ! call eigfun (nfi,iprint,tbuff,trhor,trhow,& ! &c0,irb,eigrb,bec,rhovan,rhor,rhog,rhos,& ! &2,230,276) ! !scale,instate,finstate ! call mpi_finalize(i) ! stop !###################################################################### ! !======================================================================= ! ! verlet algorithm ! ! loop which updates electronic degrees of freedom ! cm=c(t+dt) is obtained from cm=c(t-dt) and c0=c(t) ! the electron mass rises with g**2 ! !======================================================================= ! call newd(rhor,tfor,irb,eigrb,rhovan,deeq,fion) call tictac(7,0) call prefor(eigr,betae) ! !==== set friction ==== ! verl1=2./(1.+frice) verl2=1.-verl1 verl3=1./(1.+frice) ! !==== start loop ==== ! do i=1,n,2 call dforce & & (bec,deeq,betae,i,c0(1,i),c0(1,i+1),c2,c3,rhos,tbuff) if(evalue.ne.0.d0) then call dforceb & & (c0, i, betae, ipolp, bec ,ctabin(1,1,ipolp), gqq, gqqm, qmat, deeq, df) do ig=1,ngw c2(ig)=c2(ig)+evalue*df(ig) enddo call dforceb & &(c0, i+1, betae, ipolp, bec ,ctabin(1,1,ipolp), gqq, gqqm, qmat, deeq, df) do ig=1,ngw c3(ig)=c3(ig)+evalue*df(ig) enddo endif if (tnosee) then do j=1,ngw cm(j, i)=cm(j, i) & & +2.*fccc*(c0(j, i)-cm(j, i) & & +0.5*dt2bye*ema0bg(j)*c2(j)) cm(j,i+1)=cm(j,i+1) & & +2.*fccc*(c0(j,i+1)-cm(j,i+1) & & +0.5*dt2bye*ema0bg(j)*c3(j)) end do else do j=1,ngw cm(j, i) = verl1*c0(j, i) & & + verl2*cm(j, i) & & + verl3*dt2bye*ema0bg(j)*c2(j) cm(j,i+1) = verl1*c0(j,i+1) & & + verl2*cm(j,i+1) & & + verl3*dt2bye*ema0bg(j)*c3(j) end do endif if (ng0.eq.2) then cm(1, i)=cmplx(real(cm(1, i)),0.0) cm(1,i+1)=cmplx(real(cm(1,i+1)),0.0) end if end do ccc=fccc*dt2bye ! !==== end of loop which updates electronic degrees of freedom ! call tictac(7,1) ! ! buffer for wavefunctions is unit 21 ! if(tbuff) rewind 21 ! !----------------------------------------------------------------------- ! contribution to fion due to lambda !----------------------------------------------------------------------- ! ! nlfq needs deeq bec ! if(tfor)then ! call nlfq(c0,deeq,eigr,bec,becdr,fion) ! do j=1,n do i=1,n ! ! interpolate new lambda at (t+dt) from lambda(t) and lambda(t-dt): ! lambdap(i,j)=2.*lambda(i,j)-lambdam(i,j) lambdam(i,j)=lambda(i,j) lambda (i,j)=lambdap(i,j) end do end do endif ! ! calphi calculates phi ! the electron mass rises with g**2 ! call calphi(c0,ema0bg,bec,betae,phi) ! call tictac(10,0) ! ! begin try and error loop for constraints (only one step!) ! nlfl needs: lambda bec becdr ! if (tfor) call nlfl(bec,becdr,lambda,fion) ! bforceion adds the force term due to electronic berry phase ! only in US-case !if(tfor) call noforce(fion, ipolp)!ATTENZIONE if(tfor .and. (evalue .ne. 0.d0)) then call bforceion(fion,tfor,ipolp, qmat,bec,becdr,gqq,evalue) endif !------------------------------------------------ ! the center of mass along ipol is kept steady !------------------------------------------------ !if(tfor) call noforce(fion, ipolp) ! !======================================================================= ! ! verlet algorithm ! ! loop which updates ionic degrees of freedom ! taup=pos(t+dt) is obtained from taum=pos(t-dt) and tau0=pos(t) ! ! guessed displacement of ions !======================================================================= ! ! if(tfor) then ! !==== set friction ==== ! verl1=2./(1.+fricp) verl2=1.-verl1 verl3=dt2/(1.+fricp) ! if (tnosep) then do is=1,nsp do ia=1,na(is) do i=1,3 fionm(i,ia,is)=fion(i,ia,is) - & & vnhp*vel(i,ia,is)*pmass(is) taup(i,ia,is)=-taum(i,ia,is)+2.*tau0(i,ia,is)+ & & iforce(i,ia,is)*dt2*fionm(i,ia,is)/pmass(is) end do end do end do else do is=1,nsp do ia=1,na(is) do i=1,3 taup(i,ia,is) = verl1*tau0(i,ia,is) & & + verl2*taum(i,ia,is) & & + verl3/pmass(is)*fion(i,ia,is)*iforce(i,ia,is) end do end do end do endif endif ! !----------------------------------------------------------------------- ! initialization with guessed positions of ions !----------------------------------------------------------------------- ! if (tfor) then ! ! phfac calculates eigr ! call phfac(taup,ei1,ei2,ei3,eigr) ! ! prefor calculates betae ! call prefor(eigr,betae) end if ! !----------------------------------------------------------------------- ! imposing the orthogonality !----------------------------------------------------------------------- ! call tictac(11,0) if (tortho) then call ortho (eigr,cm,phi,lambda, & & bigr,iter,ccc,eps,maxit,delt,bephi,becp) else call graham(betae,bec,cm) end if call tictac(11,1) ! !----------------------------------------------------------------------- ! correction to displacement of ions !----------------------------------------------------------------------- ! if(tfor) then call cofmass(taup,cdmp) !center of mass kept steady, this should be just a numerical correction do is=1,nsp do ia=1,na(is) do i=1,3 ! comment the following line for ! NO CORRECTION TO DISPLACEMENT OF IONS, SELECTIVE POLARIZATION CALCULATION ! F.G. 25-2-03 !AAA taup(i,ia,is)=taup(i,ia,is)-cdmp(i)+cdm0(i) enddo enddo enddo call cofmass(taup,cdmp) ! write(6,*) 'Center of Mass: ', cdmp endif if(iprsta.ge.3) then nnn=min(12,n) do i=1,nnn write(6,'(a,12f18.15)') ' main lambda = ',(lambda(i,j),j=1,nnn) end do write(6,*) endif ! if (tortho) call updatc(ccc,lambda,phi,bephi,becp,bec,cm) call calbec (nvb+1,nsp,eigr,cm,bec) ! call tictac(10,1) ! if(iprsta.ge.3) call dotcsc(eigr,cm) !------------------------------------------------------------------ ! ionic center of mass kept steady on direction ipol !------------------------------------------------------------------ ! call cofcharge(taup,cdzp)!ATTENZIONE ! if(tfor) then ! call cofcharge(taup,cdzp) ! !call cofcharge(taum,cdzm) ! do is=1,nsp ! do ia=1,na(is) ! do i=1,3 ! taup(i,ia,is)=taup(i,ia,is)-cdzp(i)+cdz0(i) ! tau0(i,ia,is)=tau0(i,ia,is)-cdzp(i)+cdz0(i) ! enddo ! enddo ! enddo ! call keptwfc(cm,-cdzp(1)+cdz0(1),-cdzp(2)+cdz0(2),-cdzp(3)+cdz0(3)) ! call keptwfc(c0,-cdzp(1)+cdz0(1),-cdzp(2)+cdz0(2),-cdzp(3)+cdz0(3)) ! !call calbec (nvb+1,nsp,eigr,cm,bec) ! call cofmass(taup,cdmp) ! call cofmass(taum,cdmm) ! endif ! !----------------------------------------------------------------------- ! temperature monitored and controlled !----------------------------------------------------------------------- ! ekinv=0.0 ekinp=0.0 ekinpr=0.0 ! ! ionic kinetic energy ! if(tfor) then do is=1,nsp do ia=1,na(is) do i=1,3 vel(i,ia,is)=(taup(i,ia,is)-taum(i,ia,is))/twodel ekinp=ekinp+vel(i,ia,is)*vel(i,ia,is)*pmass(is) end do end do end do ekinp=0.5*ekinp ! ! ionic temperature ! call cofmass(vel,cdmvel) do is=1,nsp do ia=1,na(is) do i=1,3 velrel=vel(i,ia,is)-cdmvel(i) ekinpr=ekinpr+velrel*velrel*pmass(is) end do end do end do ekinpr=0.5*ekinpr endif tempp=ekinpr*factem/(1.5d0*nat) ! ! fake electronic kinetic energy ! ekinc0=0.0 do i=1,n do j=1,ngw speed=cm(j,i)-c0(j,i) ekinc0=ekinc0+2.*real(conjg(speed)*speed)/ema0bg(j) end do if(ng0.eq.2) then speed=cm(1,i)-c0(1,i) ekinc0=ekinc0-real(conjg(speed)*speed) end if end do ekinc0=ekinc0*emass/dt2 #ifdef PARA call reduce(1,ekinc0) #endif ekinc =0.5*(ekinc0+ekincm) ! ! updating nose-hoover friction variables ! if(tnosep)then xnhpp=2.*xnhp0-xnhpm+2.*(dt2/qnp)*(ekinpr-gkbt/2.) vnhp =(xnhpp-xnhpm)/twodel endif if(tnosee)then xnhep=2.*xnhe0-xnhem+2.*(dt2/qne)*(ekinc-ekincw) vnhe =(xnhep-xnhem)/twodel endif ! if(mod(nfi-1,iprint).eq.0 .or. tlast) & & call eigs(nspin,nx,nupdwn,iupdwn,f,lambda) ! epot=eht+epseu+exc press=0. enthal=etot+press*omega ! acc(1)=acc(1)+ekinc acc(2)=acc(2)+ekin acc(3)=acc(3)+epot acc(4)=acc(4)+etot acc(5)=acc(5)+tempp ! econs=ekinp+ekinv+enthal econt=econs+ekinc if(tnosep)then econt=econt+0.5*qnp*vnhp*vnhp+ gkbt*xnhp0 endif if(tnosee)then econt=econt+0.5*qne*vnhe*vnhe+2.*ekincw*xnhe0 endif ! if(mod(nfi-1,iprint).eq.0.or.tfirst) then write(6,1900) end if ! ! dis=0.0 ! discm=0.0 ! if(tfor) then ! call displace(nsp,na,taup,taui,pmass,dis,discm) ! endif ! tps=nfi*delt*2.4189d-5 ! write(6,1910) nfi,ekinc,int(tempp),etot,econs,econt,frice,fricp !the old one write(6,1910) nfi,ekinc,int(tempp),etot,econs,econt ! 1900 format(/2x,'nfi',4x,'ekinc',1x,'tempp',8x,'etot', & & 7x,'econs',7x,'econt',3x,'frice',3x,'fricp') ! 1910 format(i5,1x,f8.5,1x,i5,3(1x,f11.5),1x,2f8.5) !the old one 1910 format(i5,1x,f15.12,1x,i5,3(1x,f11.5)) ! if(tfor) then #ifdef PARA ! in parallel execution, only the first nodes writes if (me.eq.1) then #endif write(77,1910) nfi,ekinc,int(tempp),etot,econs,econt,frice,fricp if(mod(nfi-1,iprint).eq.0 .or. tlast) & & write(6,1956) 'Polarizzazione :', nfi, pberryel,pberryion,pberryel+pberryion !write(77,'(i6,1x,3f9.5)') nfi, pberryel,pberryion,evalue do is=1,nsp do ia=1,na(is) write(77,'(3f10.5)') tau0(1,ia,is), tau0(2,ia,is), tau0(3,ia,is) enddo enddo ! write(77,'(9f8.5)') & ! & (((tau0(i,ia,is),i=1,3),ia=1,na(is)),is=1,nsp) ! write(78,'(i6)') nfi!ATTENZIONE linea rimossa per il caso calcolo matrice dinamica do is=1,nsp do ia=1,na(is) write(78,'(3f30.20)') fion(1,ia,is), fion(2,ia,is), fion(3,ia,is) enddo enddo ! write(78,'(9f8.5)') & ! & (((fion(i,ia,is),i=1,3),ia=1,na(is)),is=1,nsp) #ifdef PARA endif #endif ! ! new variables for next step ! do is=1,nsp do ia=1,na(is) do i=1,3 taum(i,ia,is)=tau0(i,ia,is) tau0(i,ia,is)=taup(i,ia,is) end do end do end do ! if(tnosep) then xnhpm = xnhp0 xnhp0 = xnhpp endif if(tnosee) then xnhem = xnhe0 xnhe0 = xnhep endif do i=1,3 cdmm(i)=cdm0(i) cdm0(i)=cdmp(i) end do end if ! ekincm=ekinc0 ! ! cm=c(t+dt) c0=c(t) . Swap cm and c0 (c2 work space) ! do i=1,n do j=1,ngw c2(j)=cm(j,i) cm(j,i)=c0(j,i) c0(j,i)=c2(j) end do end do ! ! now: cm=c(t) c0=c(t+dt) ! ! write on file ndw each iprint ! if(mod(nfi,iprint).eq.0.and..not.tlast) then #ifdef PARA call writepfile & #else call writefile & #endif & ( nx,nax,nsp,ndw,nfi,ngw,n,c0,cm,tau0,vol0,cdm0,cdmm, & & taum,volm,acc,taui,a1,a2,a3,alat,lambda,lambdam, & & xnhe0,xnhem,vnhe,xnhp0,xnhpm,vnhp,vel,ekincm) endif ! ! ===================================================================== ! automatic adapt of friction using grease and twall ! ===================================================================== if (tfirst) then epre = etot enow = etot endif epre=enow enow=etot if (enow .gt. (epre+0.00002)) then if (frice .lt. grease) then savee = frice savep = fricp endif if (twall) then tbump = .true. frice = 1./grease fricp = 1./greasp endif else if (tbump) then tbump = .false. frice = savee fricp = savep endif endif frice = frice * grease fricp = fricp * greasp ! ===================================================== ! tfirst=.false. ! call tictac(1,1) ! if(nfi.lt.nomore) go to 1000 !=============================end of main loop of molecular dynamics==== ! anor=1.d0/dfloat(nfi) do i=1,nacc acc(i)=acc(i)*anor end do ! write(6,1920) 1920 format(//' averaged quantities :',/, & & 9x,'ekinc',10x,'ekin',10x,'epot',10x,'etot',5x,'tempp') write(6,1930) (acc(i),i=1,nacc) 1930 format(4f14.5,f10.1) ! #ifdef PARA call print_para_times #else call print_times #endif #ifdef PARA #ifdef AIX call reopenfile(ndw) #endif call writepfile & #else call writefile & #endif & ( nx,nax,nsp,ndw,nfi,ngw,n,c0,cm,tau0,vol0,cdm0,cdmm, & & taum,volm,acc,taui,a1,a2,a3,alat,lambda,lambdam, & & xnhe0,xnhem,vnhe,xnhp0,xnhpm,vnhp,vel,ekincm) write(6,*) ' disk written ',ndw ! if(iprsta.gt.1) then write(6,*) write(6,3370)' lambda n = ',n 3370 format(26x,a,i4) do i=1,n write(6,3380) (lambda(i,j),j=1,n) 3380 format(9f8.4) end do endif if(iprsta.gt.3) then do iss=1,nspin isa=0 do is=1,nsp do ia=1,na(is) isa=isa+1 write(6,*) write(6,3390) ' rhovan(isa) ',isa 3390 format(33x,a,i4) do iv=1,nh(is) write(6,3400) (rhovan(isa,iv,jv,iss),jv=1,nh(is)) 3400 format(8f9.4) end do write(6,*) write(6,3390) ' deeq(isa) ',isa do iv=1,nh(is) write(6,3400) (deeq(isa,iv,jv,iss),jv=1,nh(is)) end do end do end do end do endif ! if(tfor) then write(6,1940) do is=1,nsp do ia=1, na(is) ! write(6,1950) is,ia,(tau0(i,ia,is),i=1,3) & ! & ,(fion(i,ia,is),i=1,3) #ifdef PARA if(ipol.eq.1) then write(6,1955) is,ia,tau0(3,ia,is),tau0(1,ia,is),tau0(2,ia,is) elseif( ipol .eq.2) then write(6,1955) is,ia,tau0(2,ia,is),tau0(3,ia,is),tau0(1,ia,is) elseif( ipol .eq. 3) then write(6,1955) is,ia,tau0(1,ia,is),tau0(2,ia,is),tau0(3,ia,is) endif #else write(6,1955) is,ia,tau0(1,ia,is),tau0(2,ia,is),tau0(3,ia,is) #endif end do end do do is=1,nsp do ia=1, na(is) #ifdef PARA if(ipol.eq.1) then write(6,1955) is,ia,fion(3,ia,is),fion(1,ia,is),fion(2,ia,is) elseif( ipol .eq.2) then write(6,1955) is,ia,fion(2,ia,is),fion(3,ia,is),fion(1,ia,is) elseif( ipol .eq. 3) then write(6,1955) is,ia,fion(1,ia,is),fion(2,ia,is),fion(3,ia,is) endif #else write(6,1955) is,ia,fion(1,ia,is),fion(2,ia,is),fion(3,ia,is) #endif enddo enddo endif 1940 format(1x,' species',' # of atom',' x-coord', & & ' y-coord',' z-coord',/) 1950 format(1x,2i5,3f10.4,2x,3f10.4) write(6,1960) 1955 format(1x,2i5,3f11.5) 1956 format(a,i4,3(3x,f10.6)) 1960 format(5x,//'====================== end cpvan ', & & '======================',//) #ifdef PARA call mpi_finalize(i) #endif stop end