pfuiread.f

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

 	subroutine pfuiread(nunit,pfi,nfun,tutype,ielem,nderiv,ntabdim,
     $                      nhamtyp,descr,extra,ntab,atumin,atumax,
     $                      datui,tu,dtu,d2tu,icode)

*	General routine to read in potential tables from PFUI files.
*
*	FE 30-sep-93, 29-oct-93

	implicit none
*	   PARAMETERS
*	      maximum number of points in the table allowed:
	integer ntabmax
	parameter (ntabmax=20001)
*	     code for impossible hamitonian type code (quite arbitrary)
	integer impossible
	parameter (impossible=-9999999)
*	   ARGUMENTS
	integer nunit,nfun,ielem(3,nfun),nderiv,ntabdim,nhamtyp,
     $          ntab(nfun)
	integer icode
	logical pfi
	character descr(3,nfun)*80,tutype(nfun)*2
	double precision extra(3,nfun),atumin(nfun),atumax(nfun),
     $                   datui(nfun)
	double precision tu(ntabdim,nfun),dtu(ntabdim,nfun),
     $                   d2tu(ntabdim,nfun)
*	   LOCAL
	integer ideriv,itab,ifun,l,l1
	character descrT(3)*80,tutypeT*2,tutypeR*2
	integer ielemT(3),nderivT,ntabT,nhamtypT
	double precision extraT(3),atuminT,atumaxT,datuiT
	double precision tab(ntabmax)
	logical consumed
	character*5 first5
	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)

*	   initialize all the natb's to zero:
	do ifun=1,nfun
	   ntab(ifun) = 0
	enddo

*	   input checks
	if ((nderiv.lt.0).or.(nderiv.gt.2)) then
	   icode = 2
	   return
	endif

	nhamtyp = impossible
*	   read file until EOF
 11	continue
	   if (pfi) then
	      if (nunit.lt.0) then
	         read(  *  ,5,err=9999,end=12) first5
	      else
	         read(nunit,5,err=9999,end=12) first5
	      endif
	   else
	      if (nunit.ge.0) then
	         read(nunit,err=9999,end=12) first5
	      endif
	   endif
	   if (first5.ne.'PFUI1') then
*	         this does not appear to be a Level 1 PFUI file
*	         and we reject it (bleah!)
	      icode = 1
	      return
	   endif

*	      read header
	   if (pfi) then
	      if (nunit.lt.0) then
	         read(  *  ,2,err=9999) tutypeT,nderivT,ielemT,nhamtypT
	         read(  *  ,1,err=9999) descrT
	         read(  *  ,3,err=9999) extraT
	         read(  *  ,4,err=9999) ntabT
	         read(  *  ,3,err=9999) atuminT,atumaxT,datuiT
	      else
	         read(nunit,2,err=9999) tutypeT,nderivT,ielemT,nhamtypT
	         read(nunit,1,err=9999) descrT
	         read(nunit,3,err=9999) extraT
	         read(nunit,4,err=9999) ntabT
	         read(nunit,3,err=9999) atuminT,atumaxT,datuiT
	      endif
	   else
	      if (nunit.ge.0) then
	         read(nunit,err=9999) tutypeT,nderivT,ielemT,nhamtypT
	         read(nunit,err=9999) descrT
	         read(nunit,err=9999) extraT
	         read(nunit,err=9999) ntabT
	         read(nunit,err=9999) atuminT,atumaxT,datuiT
	      endif
	   endif
*	      nhamtyp will be that attached to the first function
*	      encountered in the file.
	   if (nhamtyp.eq.impossible) then
	      nhamtyp = nhamtypT
	   endif

	   if (ntabT.gt.ntabmax) then
*	         there is no room in the buffer to load these tables
	      icode = 5
	      return
	   endif

	   do ideriv=0,nderivT
*	         now we are ready to read the table for this function,
*	         derivative #ideriv.  But is this of any interest?
*	         Let us see.  also note that, in principle, the same 
*	         data can go loaded in more than one place, although 
*	         it is not clear why one would ever want to do so.
	      consumed = .false.
	      do ifun=1,nfun
*	            convert tutype to uppercase:
	         do l=1,len(tutype(ifun))
	            l1 = index(lower,tutype(ifun)(l:l))
	            if (l1.gt.0) then
	               tutypeR(l:l) = upper(l1:l1)
	            else
	               tutypeR(l:l) = tutype(ifun)(l:l)
	            endif
	         enddo
	         if ( (tutypeT.eq.tutypeR) .and.
     $                (ielemT(1).eq.ielem(1,ifun)) .and.
     $                (ielemT(2).eq.ielem(2,ifun)) .and.
     $                (ielemT(3).eq.ielem(3,ifun)) ) then
*	                  YEP! but is this a derivative of interest?
	            if (nderivT.lt.nderiv) then
*	                  this a serious error condition, we shall
*	                  never be able to supply all the derivatives
*	                  that the user wants.
	               icode = 2
	               return
	            endif
	            if (ntab(ifun).eq.0) then
*	                  not already filled in a previous pass ...
	               descr(1,ifun) = descrT(1)
	               descr(2,ifun) = descrT(2)
	               descr(3,ifun) = descrT(3)
	               extra(1,ifun) = extraT(1)
	               extra(2,ifun) = extraT(2)
	               extra(3,ifun) = extraT(3)
	               ntab(ifun) = ntabT
	               atumin(ifun) = atuminT
	               atumax(ifun) = atumaxT
	               datui(ifun) = datuiT
	               if (ntabT.gt.ntabdim) then
*	                     no room in arrays!
	                  icode = 4
	                  return
	               endif
	            endif
	            if (ideriv.le.nderiv) then
*	                  yes, the user definitely wants this 
*	                  derivative, go on and read it in buffer if
*	                  it is not already there (it *may* be already
*	                  there in the unlikely case that the user
*	                  wants the same function twice).
	               if (.not.consumed) then
	                  if (pfi) then
	                     if (nunit.lt.0) then
	                        read(  *  ,3,err=9999) 
     $                              (tab(itab),itab=1,ntabT)
	                     else
	                        read(nunit,3,err=9999) 
     $                              (tab(itab),itab=1,ntabT)
	                     endif
	                  else
	                     if (nunit.ge.0) then
	                        read(nunit,err=9999) 
     $                              (tab(itab),itab=1,ntabT)
	                     endif
	                  endif
	                  consumed = .true.
	               endif
*	                  and now transfer the buffer in the 
*	                  appropriate location.
	               if (ideriv.eq.0) then
	                  do itab=1,ntabT
	                     tu(itab,ifun) = tab(itab)
	                  enddo
	               else if (ideriv.eq.1) then
	                  do itab=1,ntabT
	                     dtu(itab,ifun) = tab(itab)
	                  enddo
	               else if (ideriv.eq.2) then
	                  do itab=1,ntabT
	                     d2tu(itab,ifun) = tab(itab)
	                  enddo
	               endif
	            endif
	         endif
	      enddo
	      if (.not.consumed) then
*	            we have to move the file pointer past this
*	            table that we did not need.  While we must
*	            go thru all data anyway for pfis, we can
*	            just skip the record (much faster!) for puis.
	         if (pfi) then
	            if (nunit.lt.0) then
	               read(  *  ,3,err=9999) 
     $                     (tab(itab),itab=1,ntabT)
	            else
	               read(nunit,3,err=9999) 
     $                     (tab(itab),itab=1,ntabT)
	            endif
	         else
	            if (nunit.ge.0) then
	               read(nunit,err=9999) 
	            endif
	         endif
	      endif
	   enddo
	   goto 11

 12	continue
*	   EOF reached.
*	   now it's the time to check whether we have loaded
*	   everything the user requested:
	icode = 0
	do ifun=1,nfun
*	      we just raise an error flag if something was not found.
	   if (ntab(ifun).eq.0) icode = 3
	enddo
	return

 9999	continue
*	   error in read:
	icode = 1
	return
	end