pfuiwrit.f

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

	subroutine pfuiwrit(nunit,pfi,nhamtyp,descr,tutype,ielem,extra,
     $                      nderiv,ntab,atumin,atumax,func,icode)

*	FE 27-sep-93, 27-oct-93

	implicit none
*	   PARAMETERS
*	      maximum number of points in the table allowed:
	integer ntabmax
	parameter (ntabmax=20001)
*	   ARGUMENTS
	integer nunit,nhamtyp,ielem(3),ntab,nderiv,icode
	logical pfi
	character descr(3)*80, tutype*2
	double precision extra(3),atumin,atumax
	external func
*	   LOCAL
	integer ideriv,itab,l,l1
	double precision datu,datui,atu,ff,dff,d2ff
	double precision tab(ntabmax)
	character*26 lower,upper
	save lower,upper
	data lower / 'abcdefghijklmnopqrstuvwxyz' /
	data upper / 'ABCDEFGHIJKLMNOPQRSTUVWXYZ' /

 1	format (a80)
 2	format (a2,5i5)
 3	format (3d25.17)
 4	format (i10)
 5	format (a5)

*	   input checks
	if ((nderiv.lt.0).or.(nderiv.gt.2)) then
	   icode = 2
	   return
	endif
	if (ntab.le.1) then
	   icode = 3
	   return
	endif
	if (atumin.ge.atumax) then
	   icode = 4
	   return
	endif
	if (ntab.gt.ntabmax) then
	   icode = 5
	   return
	endif

*	   convert tutype to uppercase:
	do l=1,len(tutype)
	   l1 = index(lower,tutype(l:l))
	   if (l1.gt.0) tutype(l:l) = upper(l1:l1)
	enddo

	datu  = (atumax-atumin) / (ntab-1)
	datui = 1.d0 / datu
	if (pfi) then
	   if (nunit.lt.0) then
	      write(  *  ,5,err=9999) 'PFUI1'
	      write(  *  ,2,err=9999) tutype,nderiv,ielem,nhamtyp
	      write(  *  ,1,err=9999) descr
	      write(  *  ,3,err=9999) extra
	      write(  *  ,4,err=9999) ntab
	      write(  *  ,3,err=9999) atumin,atumax,datui
	   else
	      write(nunit,5,err=9999) 'PFUI1'
	      write(nunit,2,err=9999) tutype,nderiv,ielem,nhamtyp
	      write(nunit,1,err=9999) descr
	      write(nunit,3,err=9999) extra
	      write(nunit,4,err=9999) ntab
	      write(nunit,3,err=9999) atumin,atumax,datui
	   endif
	else
	   if (nunit.ge.0) then
	      write(nunit,err=9999) 'PFUI1'
	      write(nunit,err=9999) tutype,nderiv,ielem,nhamtyp
	      write(nunit,err=9999) descr
	      write(nunit,err=9999) extra
	      write(nunit,err=9999) ntab
	      write(nunit,err=9999) atumin,atumax,datui
	   endif
	endif

	do ideriv=0,nderiv
*	      write out table for the ideriv-th derivative
	   do itab=1,ntab
	      atu = atumin + (itab-1)*datu
	      call func(atu,ff,dff,d2ff)
	      if (ideriv.eq.0) then
	         tab(itab) = ff
	      else if (ideriv.eq.1) then
	         tab(itab) = dff
	      else
	         tab(itab) = d2ff
	      endif
	   enddo
	   if (pfi) then
	      if (nunit.lt.0) then
	         write(  *  ,3,err=9999) (tab(itab),itab=1,ntab)
	      else
	         write(nunit,3,err=9999) (tab(itab),itab=1,ntab)
	      endif
	   else
	      if (nunit.ge.0) then
	         write(nunit,err=9999) (tab(itab),itab=1,ntab)
	      endif
	   endif
	enddo
	icode = 0
	return
 9999	continue
*	   error in write:
	icode = 1
	return
	end

	subroutine pfuifeed(nderiv,ntab,atumin,atumax,tu,dtu,d2tu,
     $                      icode)
*		this is a trick to use pfuiwrit, which requires
*		a subroutine supplying the function, also when
*		(for some reason) the user already has the tables
*		computed and wants simply to write them onto a
*		pfuifile.  The user will then 'feed' the table
*		into this routine by calling pfuifeed, and later
*		call pfuiwrit, passing pfuifunc as the 'func'
*		argument.
	implicit none
*	   PARAMETER
*	      maximum number of points in the table allowed:
	integer ntabmax
	parameter (ntabmax=20001)
*	   ARGUMENTS
	integer nderiv,ntab,icode
	double precision atumin,atumax
	double precision tu(ntab),dtu(ntab),d2tu(ntab)
	double precision arg,ff,dff,d2ff
*	   LOCAL
	double precision tuL(ntabmax),dtuL(ntabmax),d2tuL(ntabmax)
	integer nderivL,ntabL
	double precision atuminL,atumaxL,datuL,datuiL
	save tuL,dtuL,d2tuL,nderivL,ntabL,atuminL,atumaxL,datuL,datuiL
	integer itab

*	   input checks
	if ((nderiv.lt.0).or.(nderiv.gt.2)) then
	   icode = 2
	   return
	endif
	if (ntab.le.1) then
	   icode = 3
	   return
	endif
	if (atumin.ge.atumax) then
	   icode = 4
	   return
	endif
	if (ntab.gt.ntabmax) then
	   icode = 5
	   return
	endif

	nderivL = nderiv
	ntabL   = ntab
	atuminL = atumin
	atumaxL = atumax
	datuL   = (atumax-atumin) / (ntab-1)
	datuiL  = 1.d0 / datuL
	do itab=1,ntab
	   tuL(itab) = tu(itab)
	   if (nderiv.ge.1) dtuL(itab) = dtu(itab)
	   if (nderiv.ge.2) d2tuL(itab) = d2tu(itab)
	enddo
	return

	entry pfuifunc(arg,ff,dff,d2ff)
*	   arg converted into the nearest index 'itab'.  usually, 
*	   pfuifunc is called with arg exactly corresponding to a table
*	   entry, so that there are no approximations involved in
*	   the code below.
	itab = 1 + nint((arg-atuminL)*datuiL)
	if ((itab.ge.1).and.(itab.le.ntabL)) then
	   ff = tuL(itab)
	   if (nderivL.ge.1) dff = dtuL(itab)
	   if (nderivL.ge.2) d2ff = d2tuL(itab)
	else
	   ff = 0.d0
	   if (nderivL.ge.1) dff = 0.d0
	   if (nderivL.ge.2) d2ff = 0.d0
	endif
	end