!----------------------------------------------------------------------- subroutine init1(tau0,tbuff,iforce) !----------------------------------------------------------------------- ! ! read from unit 10, initialize G-vectors and related quantities ! ! ATTENZIONE l'iforce diventa vettore per tutte le componenti cartesiane use gvec use ions use parm use parms use elct use cnst use parmb use control ! implicit none ! logical tbuff integer iforce(3,nax,nsx) real(kind=8) tau0(3,nax,nsx) ! integer ibrav, idum, ik, k, iss, i, in, is, ia integer good_fft_dimension, good_fft_order real(kind=8) gcut, gcutw, gcuts, gcutb, ecut, ecutw, dual, fsum, & & b1(3), b2(3), b3(3), celldm(6), & & b1b(3),b2b(3),b3b(3), alatb, taus(3) real(kind=8) delta_pos integer i_pos,isa ! ! rewind 10 ! =================================================== ! ==== cell dimensions and lattice vectors ===== ! =================================================== ! read(10,*) ibrav,(celldm(i),i=1,6) if(ibrav.eq.0) then do i=1,3 read(10,*) a1(i),a2(i),a3(i) enddo alat=sqrt(a1(1)**2+a1(2)**2+a1(3)**3) else alat=celldm(1) end if call latgen(ibrav,celldm,a1,a2,a3,omega) ! ! a1,a2,a3 are the crystal axis (the vectors generating the lattice) ! call recips(alat,a1,a2,a3,b1,b2,b3) write(6,*) write(6,210) 210 format(' unit vectors of full simulation cell',/, & & ' in real space:',25x,'in reciprocal space:') write(6,'(3f10.4,10x,3f10.4)') a1,b1 write(6,'(3f10.4,10x,3f10.4)') a2,b2 write(6,'(3f10.4,10x,3f10.4)') a3,b3 ! ! b1,b2,b3 are the 3 basis vectors generating the reciprocal lattice ! in 2pi/alat units ! do i=1,3 ainv(1,i)=b1(i)/alat ainv(2,i)=b2(i)/alat ainv(3,i)=b3(i)/alat end do ! ! ainv is transformation matrix from cartesian to crystal coordinates ! if r=x1*a1+x2*a2+x3*a3 => x(i)=sum_j ainv(i,j)r(j) ! Note that ainv is really the inverse of a=(a1,a2,a3) ! ! ========================================================= ! ==== species, states, electrons, occupations, odd/ev ==== ! ========================================================= read(10,*) nsp,n,nel,nspin if(nsp.gt.nsx) call error(' init1',' nsp.gt.nsx',nsp) if(nspin.lt.1.or.nspin.gt.2) & & call error(' init1',' wrong nspin',nspin) ! ! important: if n is odd then nx=n+1 and c(*,n+1)=0. ! if(mod(n,2).ne.0) then nx=n+1 else nx=n end if iupdwn(1)=1 if(nspin.eq.1) then nupdwn(1)=n else read(10,*) nupdwn(1) nupdwn(2)=n-nupdwn(1) iupdwn(2)=nupdwn(1)+1 if(nupdwn(1).lt.0.or.nupdwn(2).lt.0) call error('init1', & & 'negative occupation',nupdwn(2)) endif allocate(f(nx)) allocate(ispin(nx)) read(10,*) (f(i),i=1,n) do iss=1,nspin do in=iupdwn(iss),iupdwn(iss)-1+nupdwn(iss) ispin(in)=iss end do end do ! ! ========================================================= ! ==== energy cutoffs ==== ! ========================================================= read(10,*) ecutw,ecut,dual if (dual.ne.4.0) & & call error('init1','dual.ne.4 no longer implemented',1) ! ! ========================================================= ! ==== fft grids ==== ! ========================================================= read(10,*) nr1, nr2, nr3 read(10,*) nr1s,nr2s,nr3s if (nr1s.gt.nr1.or.nr2s.gt.nr2.or.nr3s.gt.nr3) & & call error('init1','smooth grid larger than dense grid?',1) read(10,*) nr1b,nr2b,nr3b if (nr1b.gt.nr1.or.nr2b.gt.nr2.or.nr3b.gt.nr3) & & call error('init1','box grid larger than dense grid ?',1) ! ! ========================================================= ! ==== other control parameters ! ========================================================= read(10,*) tbuff,iprsta ! write(6,*) ' init1: tbuff iprsta ', & & tbuff,iprsta if (iprsta .gt. 0) then write(6,*) ' ' write(6,*) ' iprsta = 0 minimum output ' write(6,*) ' iprsta = 1 standard ' write(6,*) ' iprsta = 2 more (ok for md) ' write(6,*) ' iprsta = 3 even more (not for md) ' write(6,*) ' iprsta = 4 even more (not for md) ' write(6,*) ' iprsta = 5 all ' write(6,*) ' ' endif fsum=0. do i=1,n fsum=fsum+f(i) end do write(6,*) ' init1: fsum = ',fsum if(abs(fsum-float(nel)).gt.0.001) & & call error(' init1',' sum f(i) .ne. nel',nel) ! do is=1,nsp read(10,*) idum,na(is),ipp(is),zv(is),rcmax(is),pmass(is) if(na(is).gt.nax) call error(' init1',' na.gt.nax',na(is)) do ia=1,na(is) if(ibrav.eq.0) then read(10,*) (taus(i),i=1,3),(iforce(i,ia,is),i=1,3) do i=1,3 tau0(i,ia,is)=a1(i)*taus(1)+a2(i)*taus(2)+a3(i)*taus(3) end do else read(10,*) (tau0(i,ia,is),i=1,3),(iforce(i,ia,is),i=1,3) end if end do end do write(6,*) 'eccomi!' ! 8658 format (i5,f10.4) 8654 format (i5,3f10.4) ! #ifdef PHONON !parte per calcolare i fononi read(10, *) delta_pos read(10, *) i_pos!indice combinato if(i_pos .ne. 0) then isa = 0 do is=1,nsp do ia=1,na(is) isa=isa+1 if(i_pos .gt. 0) then if( i_pos .eq. ((isa-1)*3 +1)) then tau0(1, ia, is) = tau0(1, ia, is) + delta_pos endif if( i_pos .eq. ((isa-1)*3 +2)) then tau0(2, ia, is) = tau0(2, ia, is) + delta_pos endif if( i_pos .eq. ((isa-1)*3 +3)) then tau0(3, ia, is) = tau0(3, ia, is) + delta_pos endif else if( (-i_pos) .eq. ((isa-1)*3 +1)) then tau0(1, ia, is) = tau0(1, ia, is) - delta_pos endif if( (-i_pos) .eq. ((isa-1)*3 +2)) then tau0(2, ia, is) = tau0(2, ia, is) - delta_pos endif if( (-i_pos) .eq. ((isa-1)*3 +3)) then tau0(3, ia, is) = tau0(3, ia, is) - delta_pos endif endif enddo enddo end if #endif ! ! background electrons ! qbac=0. do is=1,nsp qbac=qbac+na(is)*zv(is) end do qbac=qbac-nel ! ! ============================================================== ! ==== set parameters and cutoffs ==== ! ============================================================== ! tpiba=2.d0*pi/alat tpiba2=tpiba*tpiba gcutw=ecutw/tpiba/tpiba gcuts=dual*gcutw gcut =ecut/tpiba/tpiba ! ! =================================================== ! ==== initialization for fft ==== ! =================================================== ! ! dense grid ! call set_fft_grid(a1,a2,a3,alat,gcut,nr1,nr2,nr3) ! ! smooth grid ! call set_fft_grid(a1,a2,a3,alat,gcuts,nr1s,nr2s,nr3s) #ifdef PARA call set_fft_para ( b1, b2, b3, gcut, gcuts, gcutw, & & nr1, nr2, nr3, nr1s, nr2s, nr3s, nnr, & & nr1x,nr2x,nr3x,nr1sx,nr2sx,nr3sx, nnrs ) #else call set_fft_dim ( nr1, nr2, nr3, nr1s, nr2s, nr3s, nnr, & & nr1x,nr2x,nr3x,nr1sx,nr2sx,nr3sx,nnrs ) #endif ! ! box grid ! nr1b=good_fft_order(nr1b) nr2b=good_fft_order(nr2b) nr3b=good_fft_order(nr3b) nr1bx=good_fft_dimension(nr1b) nr2bx=nr2b nr3bx=nr3b ! nnrb=nr1bx*nr2bx*nr3bx ! ! ============================================================== ! ==== generate g-space ==== ! ============================================================== call ggen & & ( b1, b2, b3, nr1, nr2, nr3, nr1s, nr2s, nr3s, & & nr1x, nr2x, nr3x, nr1sx, nr2sx, nr3sx, gcut, gcuts, gcutw ) ! ! ============================================================== ! generation of little box g-vectors ! ============================================================== ! alatb=alat/nr1*nr1b tpibab=2.d0*pi/alatb gcutb=ecut/tpibab/tpibab do i=1,3 a1b(i)=a1(i)/nr1*nr1b a2b(i)=a2(i)/nr2*nr2b a3b(i)=a3(i)/nr3*nr3b end do omegab=omega/nr1*nr1b/nr2*nr2b/nr3*nr3b ! call recips(alatb,a1b,a2b,a3b,b1b,b2b,b3b) write(6,*) write(6,220) 220 format(' unit vectors of full simulation cell',/, & & ' in real space:',25x,'in reciprocal space:') write(6,'(3f10.4,10x,3f10.4)') a1,b1 write(6,'(3f10.4,10x,3f10.4)') a2,b2 write(6,'(3f10.4,10x,3f10.4)') a3,b3 do i=1,3 ainvb(1,i)=b1b(i)/alatb ainvb(2,i)=b2b(i)/alatb ainvb(3,i)=b3b(i)/alatb end do ! call ggenb ( b1b, b2b, b3b, & & nr1b, nr2b, nr3b, nr1bx, nr2bx, nr3bx, gcutb) ! ! ============================================================== ! write(6,34) ibrav,alat,omega,gcut,gcuts,gcutw,1 write(6,81) nr1, nr2, nr3, nr1x, nr2x, nr3x, & & nr1s,nr2s,nr3s,nr1sx,nr2sx,nr3sx, & & nr1b,nr2b,nr3b,nr1bx,nr2bx,nr3bx ! write(6,38) write(6,334) ecutw,dual*ecutw,ecut ! write(6,6) nel,n,qbac if(nspin.eq.1)then write(6,166) nspin write(6,74) write(6,77) (f(i),i=1,n) else write(6,167) nspin,nupdwn(1),nupdwn(2) write(6,75) write(6,77) (f(i),i=iupdwn(1),nupdwn(1)) write(6,76) write(6,77) (f(i),i=iupdwn(2),iupdwn(2)-1+nupdwn(2)) endif write(6,878) nsp do is=1,nsp write(6,33) is,na(is),int(zv(is)),pmass(is),rcmax(is) write(6,9) do ia=1,na(is),3 write(6,555) ((tau0(k,ik,is),k=1,3),ik=ia,min(ia+2,na(is)) ) end do 555 format(3(4x,3(1x,f6.2))) end do ! ! ewald coeff. are transformed into cut-off radii ! do is=1,nsp rcmax(is)=1.d0/rcmax(is)**0.5d0 pmass(is)=pmass(is)*scmass end do ! 1 format(2i5,6f10.5,/,8i5,2f6.2) 2 format(3f16.10,6x,f16.10) 3 format(2i5,6f10.6) 4 format(3f24.10) 5 format(3f10.2) 33 format(' is=',i3,/,' na=',i4,' zval=',i4, & & ' atomic mass=',f6.2,' gaussian rcmax=',f6.2) 331 format(3(3(f7.4,2f16.7),/)) 34 format(' initialization ',//, & & ' ibrav=',i3,' alat=',f7.3,' omega=',f10.4, & & /,' gcut=',f7.2,3x,' gcuts=', & & f7.2,' gcutw=',f7.2,/, & & ' k-points: nkpt=',i2,//) 81 format(' meshes:',/, & & ' dense grid: nr1 ,nr2, nr3 = ',3i4, & & ' nr1x, nr2x, nr3x = ',3i4,/, & & ' smooth grid: nr1s,nr2s,nr3s = ',3i4, & & ' nr1sx,nr2sx,nr3sx= ',3i4,/, & & ' box grid: nr1b,nr2b,nr3b = ',3i4, & & ' nr1bx,nr2bx,nr3bx= ',3i4,/) 6 format(/' # of electrons=',i5,' # of states=',i5, & & ' background charge ',f5.2,/) 38 format(' explicit calculation of xc',/) 334 format(' ecutw=',f7.1,' ryd',3x, & & ' ecuts=',f7.1,' ryd',3x,' ecut=',f7.1,' ryd') 35 format(3(1x,f6.2),6x,f6.4,14x,i5) 7 format(25f3.0) 166 format(/,' nspin=',i2) 167 format(/,' nspin=',i2,5x,' nup=',i5,5x,' ndown=',i5) 74 format(' occupation numbers:') 75 format(' occupation numbers up:') 76 format(' occupation numbers down:') 77 format(20f4.1) 878 format(/' # of atomic species',i5) 9 format(' atomic coordinates:') ! return end