pfuiman.f

一个有用的分子动力学小辅助程序

	program pfuiman

*	manipulation of potentials stored in PFUI format.
*	potentials may be converted from one type to another,
*	and rescaled (if needed) during the process.
*
*	FE 12-nov-1993, 1-dec-1993

	implicit none
*** input potential data structure
	integer nfunmax,ntabmax
	parameter (nfunmax=5,ntabmax=20001)
	integer ielemI(3,nfunmax),ntabI(nfunmax),nfunI,nderivI
	double precision tuI(ntabmax,nfunmax),dtuI(ntabmax,nfunmax),
     $                   d2tuI(ntabmax,nfunmax)
	double precision extraI(3,nfunmax),atuminI(nfunmax),
     $                   atumaxI(nfunmax),datuiI(nfunmax)
	common/potI/tuI,dtuI,d2tuI,extraI,atuminI,atumaxI,datuiI,
     $              ielemI,ntabI,nfunI,nderivI
	character*2  tutypeI(nfunmax)
	character*80 descrI(3,nfunmax)
	character*10 elementI
	common/posI/tutypeI,descrI,elementI
	save /potI/,/posI/
***
	integer leneff
	external leneff
	double precision alatI,alatO,enerI,enerO,scalen,scaene,amass
	integer li,lo,icode,nhamI,nhamO,ielementO
	logical pfiI,pfiO
	character*80 fileI,fileO
	character*10 elemreq
 1	format (1x,a)
 2	format (1x,a,$)

	print '(/,1x,a,12x,a,/)','pfuiman  1-dec-93',
     $        'General purpose utility program for PFUI files'
 9991	continue
	print 2,'Input file name (must end with .pfi or .pui) ? '
	read(*,'(a)',end=9999) fileI
	li = leneff(fileI)
	if (li.lt.5) goto 9991
	if (fileI(li-3:li).eq.'.pfi') then
	   pfiI = .true.
	else if (fileI(li-3:li).eq.'.pui') then
	   pfiI = .false.
	else
	   goto 9991
	endif

	call pfuiin(fileI(1:li),pfiI,nhamI,icode)
	if (icode.ne.0) then
	   print '(1x,a,i5,a)',
     $           'pfuiin: error code',icode,' from pfuiread.'
	   goto 9991
	endif

 3300	continue
	print '(1x,3a,$)','Output element (<RETURN>=',
     $        elementI(1:leneff(elementI)),') ? '
	read(*,'(a)',end=9999) elemreq
	if (elemreq.ne.' ') then
	   ielementO = 0
	   call mendelev(ielementO,elemreq,amass)
	   if (ielementO.eq.0) then
	      print '(1x,3a)','What the heck is ',
     $              elemreq(1:leneff(elemreq)),'?!? I do not know it.'
	      goto 3300
	   endif
*	      this fixes the case:
	   elemreq = ' '
	   call mendelev(ielementO,elemreq,amass)
	   print '(1x,3a,i2,a,f9.4)',
     $           'Oh, I know ',elemreq(1:leneff(elemreq)),
     $           ': it has Z=',ielementO,', atomic mass',amass
	   print '(1x,2a)','Some numbers are needed to compute ',
     $                     'the length and energy scaling ratios:'
	   print '(1x,2a,$)',elementI(1:leneff(elementI)),
     $           ' reference length ? '
	   read(*,*,end=9999) alatI
	   print '(1x,2a,$)',elemreq(1:leneff(elemreq)),
     $           ' reference length ? '
	   read(*,*,end=9999) alatO
	   if ((alatO.eq.alatI).or.(alatI.eq.0.d0)) then
	      scalen = 0.d0
	   else
	      scalen = alatO/alatI
	   endif
	   print '(1x,2a,$)',elementI(1:leneff(elementI)),
     $           ' reference energy ? '
	   read(*,*,end=9999) enerI
	   print '(1x,2a,$)',elemreq(1:leneff(elemreq)),
     $           ' reference energy ? '
	   read(*,*,end=9999) enerO
	   if ((enerO.eq.enerI).or.(enerI.eq.0.d0)) then
	      scaene = 0.d0
	   else
	      scaene = enerO/enerI
	   endif
	   call rescale(ielementO,scalen,scaene)
	endif
 9993	continue
	print*,'The following conversions are supported at present:'
	print*,'                   11 -> 11, 12, 21'
	print*,'                   12 -> 11, 12'
	print*,'                   21 -> 21'
	print*,'                   22 -> 22'
	print 2,
     $       'Type for output potential [11 12 21 22] ? '
	read(*,*,end=9999) nhamO
	call convert(nhamI,nhamO,icode)
	if (icode.ne.0) goto 9993
 9992	continue
	print 2,'Output file name (must end with .pfi or .pui) ? '
	read(*,'(a)',end=9999) fileO
	lo = leneff(fileO)
	if (lo.lt.5) goto 9992
	if (fileO(lo-3:lo).eq.'.pfi') then
	   pfiO = .true.
	else if (fileO(lo-3:lo).eq.'.pui') then
	   pfiO = .false.
	else
	   goto 9992
	endif

	call pfuiout(fileO(1:lo),pfiO,nhamO,icode)
	if (icode.ne.0) then
	   print '(1x,a,i5,a)',
     $           'pfuiout: error code',icode,' from pfuiwrit.'
	endif
	print*,'Done.'
 9999	continue
	end

	subroutine rescale(ielementO,scalen,scaene)
*		rescale length and energy scale of input potential
	implicit none
	integer ielementO
	double precision scalen,scaene
*** input potential data structure
	integer nfunmax,ntabmax
	parameter (nfunmax=5,ntabmax=20001)
	integer ielemI(3,nfunmax),ntabI(nfunmax),nfunI,nderivI
	double precision tuI(ntabmax,nfunmax),dtuI(ntabmax,nfunmax),
     $                   d2tuI(ntabmax,nfunmax)
	double precision extraI(3,nfunmax),atuminI(nfunmax),
     $                   atumaxI(nfunmax),datuiI(nfunmax)
	common/potI/tuI,dtuI,d2tuI,extraI,atuminI,atumaxI,datuiI,
     $              ielemI,ntabI,nfunI,nderivI
	character*2  tutypeI(nfunmax)
	character*80 descrI(3,nfunmax)
	character*10 elementI
	common/posI/tutypeI,descrI,elementI
	save /potI/,/posI/
***
	integer j,itab,ifun,itry,itype(nfunmax)
	double precision factor,datu
	logical dolen,doene,quit
*	   the following tables contain, for each possible function,
*	   its scaling exponent with respect to length and energy
*	   rescaling
	integer ntypes
	parameter (ntypes=7)
	character*2 tutype(ntypes)
	integer nsclen(ntypes),nscene(ntypes)
	save tutype,nsclen,nscene
	data tutype /'V2','UU','ZT','PI','RH','FN','GN'/
	data nsclen /  1 ,  0 ,  0 ,  1 ,  1 ,  1 ,  0 /
	data nscene /  1 ,  1 ,  0 ,  1 ,  0 ,  0 ,  0 /

	dolen = (scalen.ne.0.d0)
	doene = (scaene.ne.0.d0)	
	if ((.not.dolen).and.(.not.doene)) return

*	   check we know the scaling behavior of all the functions:
	quit = .false.
	do ifun=1,nfunI
	   itype(ifun) = 0
	   do itry=1,ntypes
	      if (tutypeI(ifun).eq.tutype(itry)) then
	         itype(ifun) = itry
	         goto 16
	      endif
	   enddo
	   print '(1x,2a)',
     $           'rescale: do not know how to rescale ',tutype(ifun)
	   quit = .true.
 16	   continue
	enddo

	if (quit) then
	   print '(1x,a)','No rescaling at all was done.'
	   return
	endif

	do ifun=1,nfunI
	   do j=1,3
	      if (ielemI(j,ifun).ne.0) ielemI(j,ifun) = ielementO
	   enddo
	enddo

	if (dolen) then
*	      length scaling
	   do ifun=1,nfunI
	      if (nsclen(itype(ifun)).ne.0) then
	         factor = scalen**nsclen(itype(ifun))
	         atuminI(ifun) = atuminI(ifun)*factor
	         atumaxI(ifun) = atumaxI(ifun)*factor
	         datu = ( atumaxI(ifun) - atuminI(ifun) ) /
     $                         ( ntabI(ifun) - 1 )
	         datuiI(ifun) = 1.d0/datu
	         do itab=1,ntabI(ifun)
	            if (nderivI.ge.1) 
     $                 dtuI(itab,ifun) = dtuI(itab,ifun)/factor
	            if (nderivI.ge.2) 
     $                 d2tuI(itab,ifun) = d2tuI(itab,ifun)/(factor**2)
	         enddo
	         print '(1x,3a,f18.14)',
     $                 'Argument scale  of ',tutypeI(ifun),
     $                 ' rescaled by',factor
	      endif
	   enddo
	endif

	if (doene) then
*	      energy scaling
	   do ifun=1,nfunI
	      if (nscene(itype(ifun)).ne.0) then
	         factor = scaene**nscene(itype(ifun))
	         do itab=1,ntabI(ifun)
	            tuI(itab,ifun) = tuI(itab,ifun)*factor
	            if (nderivI.ge.1) 
     $                 dtuI(itab,ifun) = dtuI(itab,ifun)*factor
	            if (nderivI.ge.2) 
     $                 d2tuI(itab,ifun) = d2tuI(itab,ifun)*factor
	         enddo
	         print '(1x,3a,f18.14)',
     $                 'Function values of ',tutypeI(ifun),
     $                 ' rescaled by',factor
	      endif
	   enddo
	endif
	end

	subroutine convert(nhamI,nhamO,icode)
*		convert hamiltonian type nhamI into hamiltonian
*		type nhamO (just the driving logic here).
	implicit none
	integer nhamI,nhamO,icode
	logical unsupp
	unsupp = .false.
	if (nhamO.eq.nhamI) then
	   call copy
	else
	   if ((nhamI.eq.11).or.(nhamI.eq.1)) then
*	         input is glue, regular coordination
	      if (nhamO.eq.12) then
	         call conv1112
	      else if (nhamO.eq.21) then
	         call conv1121
	      else
	         unsupp = .true.
	      endif
	   else if (nhamI.eq.12) then
*	         input is glue, angular coordination
	      if (nhamO.eq.11) then
	         call conv1211
	      else
	         unsupp = .true.
	      endif
	   else
              unsupp = .true.
	   endif
	endif
	if (unsupp) then
	   print '(1x,a,i4,a,i4,a)',
     $           'Sorry, conversion from',nhamI,' to',nhamO,
     $           ' is not supported.'
	   icode = 1
	else
	   icode = 0
	endif
	end

	subroutine copy
*		copy input potential into output potential with
*		no changes whatsoever.
	implicit none
*** input potential data structure
	integer nfunmax,ntabmax
	parameter (nfunmax=5,ntabmax=20001)
	integer ielemI(3,nfunmax),ntabI(nfunmax),nfunI,nderivI
	double precision tuI(ntabmax,nfunmax),dtuI(ntabmax,nfunmax),
     $                   d2tuI(ntabmax,nfunmax)
	double precision extraI(3,nfunmax),atuminI(nfunmax),
     $                   atumaxI(nfunmax),datuiI(nfunmax)
	common/potI/tuI,dtuI,d2tuI,extraI,atuminI,atumaxI,datuiI,
     $              ielemI,ntabI,nfunI,nderivI
	character*2  tutypeI(nfunmax)
	character*80 descrI(3,nfunmax)
	character*10 elementI
	common/posI/tutypeI,descrI,elementI
	save /potI/,/posI/
***
*** output potential data structure
C	integer nfunmax,ntabmax
C	parameter (nfunmax=5,ntabmax=20001)
	integer ielemO(3,nfunmax),ntabO(nfunmax),nfunO,nderivO
	double precision tuO(ntabmax,nfunmax),dtuO(ntabmax,nfunmax),
     $                   d2tuO(ntabmax,nfunmax)
	double precision extraO(3,nfunmax),atuminO(nfunmax),
     $                   atumaxO(nfunmax),datuiO(nfunmax)
	common/potO/tuO,dtuO,d2tuO,extraO,atuminO,atumaxO,datuiO,
     $              ielemO,ntabO,nfunO,nderivO
	character*2  tutypeO(nfunmax)
	character*80 descrO(3,nfunmax)
	character*10 elementO
	common/posO/tutypeO,descrO,elementO
	save /potO/,/posO/
***
*	   LOCAL
	integer i,ifun,itab

	nfunO   = nfunI
	nderivO = nderivI
	do ifun=1,nfunO
	   tutypeO(ifun) = tutypeI(ifun)
	   do i=1,3
	      ielemO(i,ifun) = ielemI(i,ifun)
	      descrO(i,ifun) = descrI(i,ifun)
	      extraO(i,ifun) = extraI(i,ifun)
	   enddo
	   ntabO(ifun) = ntabI(ifun)
	   atuminO(ifun) = atuminI(ifun)
	   atumaxO(ifun) = atumaxI(ifun)
	   datuiO(ifun)  = datuiI(ifun)
	   do itab=1,ntabO(ifun)
	      tuO(itab,ifun) = tuI(itab,ifun)
	      if (nderivO.ge.1) dtuO(itab,ifun) = dtuI(itab,ifun)
	      if (nderivO.ge.2) d2tuO(itab,ifun) = d2tuI(itab,ifun)
	   enddo
	enddo
	end

	subroutine conv1112
*		performs conversion from 11 to 12.
*		this amounts to defining fn=sqrt(rh), and gn=0
*		everywhere except close to +1 (theta=0), where gn=1.
	implicit none
C*	   gaussian width for GN
C	double precision sigma
C	parameter (sigma=0.01d0)
*** input potential data structure
	integer nfunmax,ntabmax
	parameter (nfunmax=5,ntabmax=20001)
	integer ielemI(3,nfunmax),ntabI(nfunmax),nfunI,nderivI
	double precision tuI(ntabmax,nfunmax),dtuI(ntabmax,nfunmax),
     $                   d2tuI(ntabmax,nfunmax)
	double precision extraI(3,nfunmax),atuminI(nfunmax),
     $                   atumaxI(nfunmax),datuiI(nfunmax)
	common/potI/tuI,dtuI,d2tuI,extraI,atuminI,atumaxI,datuiI,
     $              ielemI,ntabI,nfunI,nderivI
	character*2  tutypeI(nfunmax)
	character*80 descrI(3,nfunmax)
	character*10 elementI
	common/posI/tutypeI,descrI,elementI
	save /potI/,/posI/
***
*** output potential data structure
C	integer nfunmax,ntabmax
C	parameter (nfunmax=5,ntabmax=20001)
	integer ielemO(3,nfunmax),ntabO(nfunmax),nfunO,nderivO
	double precision tuO(ntabmax,nfunmax),dtuO(ntabmax,nfunmax),
     $                   d2tuO(ntabmax,nfunmax)
	double precision extraO(3,nfunmax),atuminO(nfunmax),
     $                   atumaxO(nfunmax),datuiO(nfunmax)
	common/potO/tuO,dtuO,d2tuO,extraO,atuminO,atumaxO,datuiO,
     $              ielemO,ntabO,nfunO,nderivO
	character*2  tutypeO(nfunmax)
	character*80 descrO(3,nfunmax)
	character*10 elementO
	common/posO/tutypeO,descrO,elementO
	save /potO/,/posO/
***
*	   LOCAL
	integer i,ifun,itab
	double precision datu,arg,x2,ff,dff,d2ff

	nfunO   = 4
	nderivO = nderivI
	do ifun=1,nfunI
	   tutypeO(ifun) = tutypeI(ifun)
	   do i=1,3
	      ielemO(i,ifun) = ielemI(i,ifun)
	      descrO(i,ifun) = descrI(i,ifun)
	      extraO(i,ifun) = extraI(i,ifun)
	   enddo
	   ntabO(ifun) = ntabI(ifun)
	   atuminO(ifun) = atuminI(ifun)
	   atumaxO(ifun) = atumaxI(ifun)
	   datuiO(ifun)  = datuiI(ifun)
	   do itab=1,ntabO(ifun)
	      tuO(itab,ifun) = tuI(itab,ifun)
	      if (nderivO.ge.1) dtuO(itab,ifun) = dtuI(itab,ifun)
	      if (nderivO.ge.2) d2tuO(itab,ifun) = d2tuI(itab,ifun)
	   enddo
*	      BUT:
	   if (tutypeI(ifun).eq.'RH') then
	      tutypeO(ifun) = 'FN'
	      descrO(3,ifun) = 'Converted by Pfuiman, FN=sqrt(RH)'
	      do itab=1,ntabO(ifun)
	         arg = tuI(itab,ifun)
	         if (arg.le.0.d0) then
	            if (arg.lt.0.d0) then
	               print '(1x,a,d12.4,a,i5,a)',
     $                       'conv1112: WARNING: RH =',
     $                       arg,' < 0 at itab =',itab,', set to 0'
	            endif
	            tuO(itab,ifun) = 0.d0
	            if (nderivO.ge.1) dtuO(itab,ifun) = 0.d0
	            if (nderivO.ge.2) d2tuO(itab,ifun) = 0.d0
	         else
	            tuO(itab,ifun) = sqrt(arg)
	            if (nderivO.ge.1) dtuO(itab,ifun) = 
     $                           0.5d0*dtuI(itab,ifun)/tuO(itab,ifun)
	            if (nderivO.ge.2) d2tuO(itab,ifun) =
     $              ( 0.5d0*d2tuI(itab,ifun) - (dtuO(itab,ifun)**2) )/
     $                                tuO(itab,ifun)
	         endif
	      enddo
	   endif
	enddo
	tutypeO(4) = 'GN'
	ielemO(1,4) = ielemI(1,1)
	ielemO(2,4) = 0
	ielemO(3,4) = 0
	extraO(1,4) = 1.d0
	extraO(2,4) = 0.d0
	extraO(3,4) = 0.d0
	do i=1,3
	   descrO(i,4) = descrI(i,1)
	enddo
	descrO(3,4) = 'Null GN(cos(theta)) generated by Pfuiman.'
	ntabO(4) = ntabI(1)
	atuminO(4) = -1.d0
	atumaxO(4) = +1.d0
	datu       = (atumaxO(4) - atuminO(4)) / (ntabO(4)-1)
	datuiO(4)  = 1.d0/datu
C*	   GN is a gaussian with width sigma
	do itab=1,ntabO(4)
C	   arg = atuminO(4) + (itab-1)*datu
C	   x2 = 0.5d0*( (arg-1.d0)/sigma )**2
C	   if (x2.gt.50.d0) then
	      ff = 0.d0
	      dff = 0.d0
	      d2ff = 0.d0
C	   else
C	      ff = exp(-x2)
C	      dff = -ff*(arg-1.d0)/sigma**2
C	      d2ff = -(dff*(arg-1.d0)+ff)/sigma**2
C	   endif
	   tuO(itab,4) = ff
	   if (nderivO.ge.1) dtuO(itab,4) = dff
	   if (nderivO.ge.2) d2tuO(itab,4) = d2ff
	enddo
	end

	subroutine conv1121
*		performs conversion from 11 to 21.
*		this amounts to defining zt=sqrt(uu/u0), and 
*		psi(r)=-psi0*exp(-r^2/2sigma^2), with psi0/2=u0.
*		here, u0 is the value of uu for a perfect bulk
*		crystal.  for the same crystal, therefore, zt
*		will be normalized to 1.
*		note that we assume uu<0 (and u0<0), psi0<0, zt>0.