pfuiman.f

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

	implicit none
*	   gaussian width for PI, in angstroms
	double precision sigma
	parameter (sigma=0.1d0)
*** 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,ifunO,iv2,iuu,irh,itab,mmm
	double precision datu,arg,x2,ff,dff,d2ff,coord0,uu0,psi0,
     $                   ccc,ddd

	nfunO   = 4
	nderivO = nderivI
*	   we start with a copy.
	do ifun=1,nfunI
	   if (tutypeI(ifun).eq.'V2') then
	      iv2 = ifun
	      ifunO = 1
	   else if (tutypeI(ifun).eq.'UU') then
	      iuu = ifun
	      ifunO = 3
	   else if (tutypeI(ifun).eq.'RH') then
	      irh = ifun
	      ifunO = 4
	   else
*	         this should never happen
	      print*,'WARNING!!! tutypeI = ',tutypeI(ifun)
	      ifunO = 2
	   endif
	   tutypeO(ifunO) = tutypeI(ifun)
	   do i=1,3
	      ielemO(i,ifunO) = ielemI(i,ifun)
	      descrO(i,ifunO) = descrI(i,ifun)
	      extraO(i,ifunO) = extraI(i,ifun)
	   enddo
	   ntabO(ifunO) = ntabI(ifun)
	   atuminO(ifunO) = atuminI(ifun)
	   atumaxO(ifunO) = atumaxI(ifun)
	   datuiO(ifunO)  = datuiI(ifun)
	   do itab=1,ntabO(ifunO)
	      tuO(itab,ifunO) = tuI(itab,ifun)
	      if (nderivO.ge.1) dtuO(itab,ifunO) = dtuI(itab,ifun)
	      if (nderivO.ge.2) d2tuO(itab,ifunO) = d2tuI(itab,ifun)
	   enddo
	enddo

	print '(1x,2a,$)',
     $        'Enter reference value for coordination ',
     $        '(usually n=1): n = '
	read*,coord0
*	   interpolate in UU table to get value of glue function
*	   at this coordination:
	if (coord0.le.atuminO(3)) then
	   uu0 = 0.d0
	   print '(1x,a,g17.10)','WARNING: table starts at',atuminO(3)
	else if (coord0.lt.atumaxO(3)) then
	   ccc = (coord0 - atuminO(3))*datuiO(3) + 1.d0
	   mmm = ccc
	   if (mmm.lt.1) mmm = 1
	   if (mmm.ge.ntabO(3)) mmm = ntabO(3) - 1
	   ddd = ccc - mmm
	   uu0 = tuO(mmm+1,3)*ddd + tuO(mmm,3)*(1.d0-ddd)
	else
	   ccc = (coord0 - atuminO(3))*datuiO(3) + 1.d0
	   mmm = ntabO(3) - 1
	   ddd = ccc - mmm
	   uu0 = tuO(mmm+1,3)*ddd + tuO(mmm,3)*(1.d0-ddd)
	   print '(1x,a,g17.10)','WARNING: table ends at',atumaxO(3)
	endif
	print '(1x,a,g17.10)','Here, U = ',uu0

*	   now turn UU into ZT by taking its square root
	ifunO = 3
	tutypeO(ifunO) = 'ZT'
	descrO(3,ifunO) = 'Converted by Pfuiman, ZT=sqrt(UU/U0)'
	do itab=1,ntabO(ifunO)
	   arg = tuI(itab,iuu)/uu0
	   if (arg.le.0.d0) then
	      if (arg.lt.0.d0) then
	               print '(1x,a,d12.4,a,i5,a)',
     $                       'conv1121: WARNING: UU/U0 =',
     $                       arg,' < 0 at itab =',itab,', set to 0'
	      endif
	      tuO(itab,ifunO) = 0.d0
	      if (nderivO.ge.1) dtuO(itab,ifunO) = 0.d0
	      if (nderivO.ge.2) d2tuO(itab,ifunO) = 0.d0
	   else
	      tuO(itab,ifunO) = sqrt(arg)
	      if (nderivO.ge.1) dtuO(itab,ifunO) = 
     $                 0.5d0*(dtuI(itab,iuu)/uu0)/tuO(itab,ifunO)
	      if (nderivO.ge.2) d2tuO(itab,ifunO) =
     $        ( 0.5d0*(d2tuI(itab,iuu)/uu0) - (dtuO(itab,ifunO)**2) )/
     $                              tuO(itab,ifunO)
	   endif
	enddo

*	   so far so good, now we define the psi function as a
*	   gaussian.
*	   at r=0 it will be (that is a negative number!):
	psi0 = uu0 + uu0

	ifunO = 2
	tutypeO(ifunO) = 'PI'
	ielemO(1,ifunO) = ielemI(1,iv2)
	ielemO(2,ifunO) = ielemI(2,iv2)
	ielemO(3,ifunO) = 0
	extraO(1,ifunO) = psi0
	extraO(2,ifunO) = 0.d0
	extraO(3,ifunO) = 0.d0
	do i=1,3
	   descrO(i,ifunO) = descrI(i,1)
	enddo
	descrO(3,ifunO) = 'Gaussian Psi generated by Pfuiman.'
*	   psi *must* start from zero:
	atuminO(ifunO) =  0.d0
*	   range of psi made same as range of rho:
	atumaxO(ifunO) = atumaxI(irh)
*	   we try to retain the same table resolution of rho
	ntabO(ifunO) = nint( 1 + (ntabI(irh)-1)*
     $      (atumaxO(ifunO)-atuminO(ifunO)) /
     $      (atumaxI(irh)  -atuminI(irh)  )   )
	ntabO(ifunO) = min(ntabO(ifunO),ntabmax)
	datu       = (atumaxO(ifunO) - atuminO(ifunO)) / 
     $                       (ntabO(ifunO)-1)
	datuiO(ifunO)  = 1.d0/datu
*	   Psi is a gaussian with width sigma
	do itab=1,ntabO(ifunO)
	   arg = atuminO(ifunO) + (itab-1)*datu
	   x2 = 0.5d0*( (arg/sigma)**2 )
	   if (x2.gt.50.d0) then
	      ff = 0.d0
	      dff = 0.d0
	      d2ff = 0.d0
	   else
	      ff = psi0*exp(-x2)
	      dff = -ff*(arg-1.d0)/sigma**2
	      d2ff = -(dff*(arg-1.d0)+ff)/sigma**2
	   endif
	   tuO(itab,ifunO) = ff
	   if (nderivO.ge.1) dtuO(itab,ifunO) = dff
	   if (nderivO.ge.2) d2tuO(itab,ifunO) = d2ff
	enddo
	end

	subroutine conv1211
*		performs conversion from 12 to 11.
*		this amounts to throwing away gn, and defining
*		rh=fn**2
	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,ifo,ifun,itab,ifn
	double precision arg,gn1

	nfunO = 3
	nderivO = nderivI
	ifo = 0
	do ifun=1,nfunI
	   if (tutypeI(ifun).eq.'GN') then
*	         we just need to retain gn1 from here.
	      print*,'GN function just thrown away.'
	      gn1 = extraI(1,ifun)
	      print '(1x,a,e25.17)','gn(1) =',gn1
	   else if (tutypeI(ifun).eq.'FN') then
*	         this will become, squared, RH.  But it will appear
*	         only at the end, because gn1 is needed, and RH must
*	         canonically stay at the end in format 11 anyway.
	      ifn = ifun
	   else
*	         this is V2 or UU which is just copied:
	      ifo = ifo + 1
	      tutypeO(ifo) = tutypeI(ifun)
	      do i=1,3
	         ielemO(i,ifo) = ielemI(i,ifun)
	         descrO(i,ifo) = descrI(i,ifun)
	         extraO(i,ifo) = extraI(i,ifun)
	      enddo
	      ntabO(ifo) = ntabI(ifun)
	      atuminO(ifo) = atuminI(ifun)
	      atumaxO(ifo) = atumaxI(ifun)
	      datuiO (ifo)  = datuiI(ifun)
	      do itab=1,ntabO(ifo)
	         tuO(itab,ifo) = tuI(itab,ifun)
	         if (nderivO.ge.1) dtuO(itab,ifo) = dtuI(itab,ifun)
	         if (nderivO.ge.2) d2tuO(itab,ifo) = d2tuI(itab,ifun)
	      enddo
	   endif
	enddo

*	   at this point we compute RH and put it in third position.
*	   data for FN are indexed by ifn, previously saved.
	ifo = ifo + 1
	tutypeO(ifo) = 'RH'
	do i=1,3
	   ielemO(i,ifo) = ielemI(i,ifn)
	   descrO(i,ifo) = descrI(i,ifn)
	   extraO(i,ifo) = extraI(i,ifn)
	enddo
	descrO(3,ifo) = 'Converted by Pfuiman, RH=FN^2'
	ntabO(ifo) = ntabI(ifn)
	atuminO(ifo) = atuminI(ifn)
	atumaxO(ifo) = atumaxI(ifn)
	datuiO (ifo)  = datuiI(ifn)

	do itab=1,ntabO(ifo)
	   arg = tuI(itab,ifn)
	   tuO(itab,ifo) = gn1*arg*arg
	   if (nderivO.ge.1) dtuO(itab,ifo) =
     $         (gn1+gn1)* arg*dtuI(itab,ifn)
	   if (nderivO.ge.2) d2tuO(itab,ifo) =
     $         (gn1+gn1)*( dtuI(itab,ifn)**2 + arg*d2tuI(itab,ifn) )
	enddo

	if (ifo.ne.nfunO) then
	   print*,'conv1211: WARNING: wrong # of functions:',ifo,
     $            ', not',nfunO
	endif
	end

	subroutine pfuiin(fileI,pfiI,nhamI,icode)
*		read in PFUI file
	implicit none
	character*(*) fileI
	logical pfiI
	integer nhamI,icode
*** 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
	integer i,ifun,le,ll
	double precision amass
	character*40 line

	if (pfiI) then 
	   open(1,file=fileI,form='formatted',status='old',err=9999)
	else
	   open(1,file=fileI,form='unformatted',status='old',err=9999)
	endif
	rewind 1

	call pfuiinfo(1,pfiI,nfunmax,
     $                nfunI,nhamI,descrI,tutypeI,ielemI,extraI,
     $                nderivI,ntabI,atuminI,atumaxI,icode)
	if (icode.ne.0) then
	   print '(1x,a,i5,a)',
     $            'pfuiin: icode =',icode,' from pfuiinfo'
	   return
	endif

*	   printouts
	print '(1x,a,i3,a,i2)',
     $        'Potential type',nhamI,
     $        ', derivative order up to',nderivI
	print '(1x,a,i3,a,/)','File contains',nfunI,' functions:'
	do ifun=1,nfunI
	   do i=1,3
	      if (ifun.gt.1) then
	         if (descrI(i,ifun).ne.descrI(i,ifun-1)) then
	            print '(1x,a)',
     $                 descrI(i,ifun)(1:leneff(descrI(i,ifun)))
	         endif
	      else
	         print '(1x,a)',
     $              descrI(i,ifun)(1:leneff(descrI(i,ifun)))
	      endif
	   enddo
	   line = tutypeI(ifun) // ' for '
	   ll = 7
	   do i=1,3
	      if (ielemI(i,ifun).ne.0) then
	         call mendelev(ielemI(i,ifun),elementI,amass)
	         le = leneff(elementI)
	         line(ll+1:ll+le) = elementI
	         ll = ll + le + 1
	      endif
	   enddo
	   print '(1x,a)',line(1:ll-1)
	   print '(1x,i6,a,g14.6,a,g14.6)',
     $           ntabI(ifun),'-points table from',
     $           atuminI(ifun),' to',atumaxI(ifun)
	   print '()'
	enddo

*	   now read it in for real:
	rewind 1
	call pfuiread(1,pfiI,nfunI,tutypeI,ielemI,nderivI,ntabmax,
     $                nhamI,descrI,extraI,ntabI,atuminI,atumaxI,datuiI,
     $                tuI,dtuI,d2tuI,icode)
	if (icode.ne.0) return
	print '(1x,a)','Load successful.'
	return
 9999	continue
	icode = 9999
	print '(1x,2a)',fileI,' not found.'
	return
	end

	subroutine pfuiout(fileO,pfiO,nhamO,icode)
*		write out PFUI file
	implicit none
	character*(*) fileO
	logical pfiO
	integer nhamO,icode
*** output potential data structure
	integer nfunmax,ntabmax
	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/
***
	integer leneff
	external leneff,pfuifunc
	integer i,ifun,le,ll
	double precision amass
	character*40 line

	if (pfiO) then 
	   open(2,file=fileO,form='formatted',status='unknown',
     $          err=9999)
	else
	   open(2,file=fileO,form='unformatted',status='unknown',
     $          err=9999)
	endif
	rewind 2

*	   printouts
	print '(1x,a)','Writing out:'
	print '(1x,a,i3,a,i2)',
     $        'Potential type',nhamO,
     $        ', derivative order up to',nderivO
	print '(1x,a,i3,a,/)','File contains',nfunO,' functions:'
	do ifun=1,nfunO
	   do i=1,3
	      if (ifun.gt.1) then
	         if (descrO(i,ifun).ne.descrO(i,ifun-1)) then
	            print '(1x,a)',
     $                 descrO(i,ifun)(1:leneff(descrO(i,ifun)))
	         endif
	      else
	         print '(1x,a)',
     $              descrO(i,ifun)(1:leneff(descrO(i,ifun)))
	      endif
	   enddo
	   line = tutypeO(ifun) // ' for '
	   ll = 7
	   do i=1,3
	      if (ielemO(i,ifun).ne.0) then
	         call mendelev(ielemO(i,ifun),elementO,amass)
	         le = leneff(elementO)
	         line(ll+1:ll+le) = elementO
	         ll = ll + le + 1
	      endif
	   enddo
	   print '(1x,a)',line(1:ll-1)
	   print '(1x,i6,a,g14.6,a,g14.6)',
     $           ntabO(ifun),'-points table from',
     $           atuminO(ifun),' to',atumaxO(ifun)
	   print '()'
	enddo

	do ifun=1,nfunO
*	      initialize pseudo-function to return existing tables:
	   call pfuifeed(nderivO,ntabO(ifun),
     $                   atuminO(ifun),atumaxO(ifun),
     $                   tuO(1,ifun),dtuO(1,ifun),d2tuO(1,ifun),icode)
	   if (icode.ne.0) then
	      print '(1x,a,i5,a)',
     $              'pfuiout: icode =',icode,' from pfuifeed.'
	      return
	   endif
	   call pfuiwrit(2,pfiO,nhamO,descrO(1,ifun),tutypeO(ifun),
     $                   ielemO(1,ifun),extraO(1,ifun),nderivO,
     $                   ntabO(ifun),atuminO(ifun),atumaxO(ifun),
     $                   pfuifunc,icode)
	   if (icode.ne.0) then
	      print '(1x,a,i5,a)',
     $              'pfuiout: icode =',icode,' from pfuiwrit.'
	      return
	   endif
	   print '(1x,a,i3)','Successful dump for function #',ifun
	enddo
	close(2)
	print '(1x,2a)',fileO,' written.'
	icode = 0
	return
 9999	continue
	icode = 9999
	print '(1x,2a)','Cannot open ',fileO
	return
	end