double precision betaT
- integer iscode,indpdb,outpdb,outmol2,iopt,nstart,nend
+ integer iscode,indpdb,outpdb,outmol2,iopt,nstart,nend,constr_dist
logical refstr,pdbref,punch_dist,print_dist,caonly,lside,
- & lprint_cart,lprint_int,from_cart,efree,from_bx,from_cx
+ & lprint_cart,lprint_int,from_cart,efree,from_bx,from_cx,
+ & with_dihed_constr
common /cntrl/ betaT,iscode,indpdb,refstr,pdbref,outpdb,outmol2,
& punch_dist,print_dist,caonly,lside,lprint_cart,lprint_int,
- & from_cart,from_bx,from_cx,efree,iopt,nstart,nend
+ & from_cart,from_bx,from_cx,efree,iopt,nstart,nend,constr_dist,
+ & with_dihed_constr
--- /dev/null
+******************************************************************
+*
+* Array dimensions for level-based conformation comparison program:
+*
+* Max. number levels of comparison
+*
+ integer maxlevel
+ PARAMETER (MAXLEVEL=3)
+*
+* Max. number of fragments at a given level of comparison
+*
+ integer maxfrag,mmaxfrag
+ PARAMETER (MAXFRAG=30,MMAXFRAG=MAXFRAG*(MAXFRAG+1)/2)
+*
+* Max. number of pieces forming a substructure to be compared
+*
+ integer maxpiece
+ PARAMETER (MAXPIECE=20)
+*
+*******************************************************************
+++ /dev/null
-INSTALL_DIR = /users/software/mpich-1.2.7p1_intel-10.1_em64_ssh
-BIN=/users/adam/ZSCOREZ/bin
-FC = ifort
-OPT = -O3 -ip -w
-OPT = -CB -g
-FFLAGS = ${OPT} -c -I. -I../src_MD_T-sccor/include_unres -I../src_MD_T-sccor -I/users/adam/MEY_MD/src_Tc-czarek -I$(INSTALL_DIR)/include
-CPPFLAGS = -DLINUX -DPGI -DSPLITELE -DPROCOR -DCRYST_BOND -DCRYST_THETA -DCRYST_SC -DMP -DMPI
-LIBS = -L$(INSTALL_DIR)/lib -lmpich ../srcWHAM-Tsccor/xdrf/libxdrf.a -g -d2 -CA -CB
-
-.c.o:
- cc -c -DLINUX -DPGI $*.c
-
-.f.o:
- ${FC} ${FFLAGS} $*.f
-
-.F.o:
- ${FC} ${FFLAGS} ${CPPFLAGS} ${FFLAGS} $*.F
-
-objects = main_clust.o arcos.o cartprint.o chainbuild.o convert.o initialize_p.o \
- matmult.o readrtns.o pinorm.o rescode.o intcor.o timing.o misc.o \
- geomout.o readpdb.o read_coords.o parmread.o probabl.o fitsq.o hc.o \
- track.o wrtclust.o srtclust.o noyes.o contact.o printmat.o \
- int_from_cart1.o energy_p_new.o icant.o proc_proc.o work_partition.o \
- setup_var.o read_ref_str.o
-
-unres_clust: $(objects)
- $(FC) ${OPT} ${objects} ${LIBS} -o ${BIN}/unres_clustMD_MPI-oldparm
-
-clean:
- /bin/rm *.o
-
-move:
- mv *.o ${OBJ}
--- /dev/null
+Makefile-MPICH-ifort
\ No newline at end of file
--- /dev/null
+INSTALL_DIR = /users/software/mpich-1.2.7p1_intel-10.1_em64_ssh
+BIN=../../../../bin/cluster
+FC = ifort
+OPT = -O3 -ip -w
+#OPT = -CB -g
+FFLAGS = ${OPT} -c -I. -Iinclude_unres -I$(INSTALL_DIR)/include
+CPPFLAGS = -DLINUX -DPGI -DSPLITELE -DPROCOR -DMP -DMPI
+LIBS = -L$(INSTALL_DIR)/lib -lmpich xdrf/libxdrf.a -g -d2 -CA -CB
+
+.c.o:
+ cc -c -DLINUX -DPGI $*.c
+
+.f.o:
+ ${FC} ${FFLAGS} $*.f
+
+.F.o:
+ ${FC} ${FFLAGS} ${CPPFLAGS} ${FFLAGS} $*.F
+
+objects = main_clust.o arcos.o cartprint.o chainbuild.o convert.o initialize_p.o \
+ matmult.o readrtns.o pinorm.o rescode.o intcor.o timing.o misc.o \
+ geomout.o readpdb.o read_coords.o parmread.o probabl.o fitsq.o hc.o \
+ track.o wrtclust.o srtclust.o noyes.o contact.o printmat.o \
+ int_from_cart1.o energy_p_new.o icant.o proc_proc.o work_partition.o \
+ setup_var.o read_ref_str.o gnmr1.o
+
+GAB: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DUNRES -DISNAN -DMP -DMPI \
+ -DSPLITELE -DLANG0 -DCRYST_BOND -DCRYST_THETA -DCRYST_SC
+GAB: $(objects) xdrf/libxdrf.a
+ $(FC) ${OPT} ${objects} ${LIBS} -o ${BIN}/unres_clustMD_MPICH-GAB.exe
+
+E0LL2Y: CPPFLAGS = -DPROCOR -DLINUX -DPGI -DUNRES -DISNAN -DMP -DMPI \
+ -DSPLITELE -DLANG0
+E0LL2Y: $(objects) xdrf/libxdrf.a
+ $(FC) ${OPT} ${objects} ${LIBS} -o ${BIN}/unres_clustMD_MPICH-E0LL2Y.exe
+
+xdrf/libxdrf.a:
+ cd xdrf && make
+
+
+clean:
+ /bin/rm -f *.o && /bin/rm -f compinfo && cd xdrf && make clean
+
+
C
implicit real*8 (a-h,o-z)
include 'DIMENSIONS'
- include 'sizesclu.dat'
include 'COMMON.SBRIDGE'
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.VAR'
include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
dimension ggg(3)
ehpb=0.0D0
-cd print *,'edis: nhpb=',nhpb,' fbr=',fbr
-cd print *,'link_start=',link_start,' link_end=',link_end
+cd write(iout,*)'edis: nhpb=',nhpb,' fbr=',fbr
+cd write(iout,*)'link_start=',link_start,' link_end=',link_end
if (link_end.eq.0) return
do i=link_start,link_end
C If ihpb(i) and jhpb(i) > NRES, this is a SC-SC distance, otherwise a
iii=ii
jjj=jj
endif
+c write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj,
+c & dhpb(i),dhpb1(i),forcon(i)
C 24/11/03 AL: SS bridges handled separately because of introducing a specific
C distance and angle dependent SS bond potential.
if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then
call ssbond_ene(iii,jjj,eij)
ehpb=ehpb+2*eij
+cd write (iout,*) "eij",eij
+ else if (ii.gt.nres .and. jj.gt.nres) then
+c Restraints from contact prediction
+ dd=dist(ii,jj)
+ if (dhpb1(i).gt.0.0d0) then
+ ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+ fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd
+c write (iout,*) "beta nmr",
+c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+ else
+ dd=dist(ii,jj)
+ rdis=dd-dhpb(i)
+C Get the force constant corresponding to this distance.
+ waga=forcon(i)
+C Calculate the contribution to energy.
+ ehpb=ehpb+waga*rdis*rdis
+c write (iout,*) "beta reg",dd,waga*rdis*rdis
+C
+C Evaluate gradient.
+C
+ fac=waga*rdis/dd
+ endif
+ do j=1,3
+ ggg(j)=fac*(c(j,jj)-c(j,ii))
+ enddo
+ do j=1,3
+ ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
+ ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
+ enddo
+ do k=1,3
+ ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
+ ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
+ enddo
else
C Calculate the distance between the two points and its difference from the
C target distance.
- dd=dist(ii,jj)
- rdis=dd-dhpb(i)
+ dd=dist(ii,jj)
+ if (dhpb1(i).gt.0.0d0) then
+ ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+ fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd
+c write (iout,*) "alph nmr",
+c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+ else
+ rdis=dd-dhpb(i)
C Get the force constant corresponding to this distance.
- waga=forcon(i)
+ waga=forcon(i)
C Calculate the contribution to energy.
- ehpb=ehpb+waga*rdis*rdis
+ ehpb=ehpb+waga*rdis*rdis
+c write (iout,*) "alpha reg",dd,waga*rdis*rdis
C
C Evaluate gradient.
C
- fac=waga*rdis/dd
+ fac=waga*rdis/dd
+ endif
cd print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd,
cd & ' waga=',waga,' fac=',fac
- do j=1,3
- ggg(j)=fac*(c(j,jj)-c(j,ii))
- enddo
+ do j=1,3
+ ggg(j)=fac*(c(j,jj)-c(j,ii))
+ enddo
cd print '(i3,3(1pe14.5))',i,(ggg(j),j=1,3)
C If this is a SC-SC distance, we need to calculate the contributions to the
C Cartesian gradient in the SC vectors (ghpbx).
- if (iii.lt.ii) then
+ if (iii.lt.ii) then
do j=1,3
ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
enddo
- endif
- do j=iii,jjj-1
+ endif
do k=1,3
- ghpbc(k,j)=ghpbc(k,j)+ggg(k)
+ ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
+ ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
enddo
- enddo
endif
enddo
ehpb=0.5D0*ehpb
do i=1,ndih_constr
itori=idih_constr(i)
phii=phi(itori)
- difi=phii-phi0(i)
+ difi=pinorm(phii-phi0(i))
if (difi.gt.drange(i)) then
difi=difi-drange(i)
edihcnstr=edihcnstr+0.25d0*ftors*difi**4
edihcnstr=edihcnstr+0.25d0*ftors*difi**4
gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
endif
-! write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii,
-! & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg)
+c write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii,
+c & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg)
enddo
-! write (iout,*) 'edihcnstr',edihcnstr
+ write (iout,*) 'edihcnstr',edihcnstr
return
end
c------------------------------------------------------------------------------
enddo
! 6/20/98 - dihedral angle constraints
edihcnstr=0.0d0
+c write (iout,*) "Dihedral angle restraint energy"
do i=1,ndih_constr
- print *,"i",i
itori=idih_constr(i)
phii=phi(itori)
- difi=phii-phi0(i)
+ difi=pinorm(phii-phi0(i))
+c write (iout,'(2i5,4f8.3,2e14.5)') i,itori,rad2deg*phii,
+c & rad2deg*difi,rad2deg*phi0(i),rad2deg*drange(i)
if (difi.gt.drange(i)) then
difi=difi-drange(i)
edihcnstr=edihcnstr+0.25d0*ftors*difi**4
gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
+c write (iout,*) 0.25d0*ftors*difi**4
else if (difi.lt.-drange(i)) then
difi=difi+drange(i)
edihcnstr=edihcnstr+0.25d0*ftors*difi**4
gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3
+c write (iout,*) 0.25d0*ftors*difi**4
endif
-! write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii,
-! & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg)
enddo
-! write (iout,*) 'edihcnstr',edihcnstr
+c write (iout,*) 'edihcnstr',edihcnstr
return
end
c----------------------------------------------------------------------------
--- /dev/null
+ double precision function gnmr1(y,ymin,ymax)
+ implicit none
+ double precision y,ymin,ymax
+ double precision wykl /4.0d0/
+ if (y.lt.ymin) then
+ gnmr1=(ymin-y)**wykl/wykl
+ else if (y.gt.ymax) then
+ gnmr1=(y-ymax)**wykl/wykl
+ else
+ gnmr1=0.0d0
+ endif
+ return
+ end
+c------------------------------------------------------------------------------
+ double precision function gnmr1prim(y,ymin,ymax)
+ implicit none
+ double precision y,ymin,ymax
+ double precision wykl /4.0d0/
+ if (y.lt.ymin) then
+ gnmr1prim=-(ymin-y)**(wykl-1)
+ else if (y.gt.ymax) then
+ gnmr1prim=(y-ymax)**(wykl-1)
+ else
+ gnmr1prim=0.0d0
+ endif
+ return
+ end
+c------------------------------------------------------------------------------
+ double precision function harmonic(y,ymax)
+ implicit none
+ double precision y,ymax
+ double precision wykl /2.0d0/
+ harmonic=(y-ymax)**wykl
+ return
+ end
+c-------------------------------------------------------------------------------
+ double precision function harmonicprim(y,ymax)
+ double precision y,ymin,ymax
+ double precision wykl /2.0d0/
+ harmonicprim=(y-ymax)*wykl
+ return
+ end
+c---------------------------------------------------------------------------------
--- /dev/null
+ integer i,j,k,l
+ double precision erij,rij,xj,yj,zj,dxi,dyi,dzi,dxj,dyj,dzj,
+ & chi1,chi2,chi12,chip1,chip2,chip12,alf1,alf2,alf12,om1,om2,om12,
+ & om1om2,chiom1,chiom2,chiom12,chipom1,chipom2,chipom12,eps1,
+ & faceps1,faceps1_inv,eps1_om12,facsig,sigsq,sigsq_om1,sigsq_om2,
+ & sigsq_om12,facp,facp_inv,facp1,eps2rt,eps2rt_om1,eps2rt_om2,
+ & eps2rt_om12,eps3rt,eom1,eom2,eom12,evdwij,eps2der,eps3der,sigder,
+ & dsci_inv,dscj_inv,gg
+ common /calc/ erij(3),rij,xj,yj,zj,dxi,dyi,dzi,dxj,dyj,dzj,
+ & chi1,chi2,chi12,chip1,chip2,chip12,alf1,alf2,alf12,om1,om2,om12,
+ & om1om2,chiom1,chiom2,chiom12,chipom1,chipom2,chipom12,eps1,
+ & faceps1,faceps1_inv,eps1_om12,facsig,sigsq,sigsq_om1,sigsq_om2,
+ & sigsq_om12,facp,facp_inv,facp1,eps2rt,eps2rt_om1,eps2rt_om2,
+ & eps2rt_om12,eps3rt,eom1,eom2,eom12,evdwij,eps2der,eps3der,sigder,
+ & dsci_inv,dscj_inv,gg(3),i,j
--- /dev/null
+C Change 12/1/95 - common block CONTACTS1 included.
+ integer ncont,ncont_ref,icont,icont_ref,num_cont,jcont
+ double precision facont,gacont
+ common /contacts/ ncont,ncont_ref,icont(2,maxcont),
+ & icont_ref(2,maxcont)
+ common /contacts1/ facont(maxconts,maxres),
+ & gacont(3,maxconts,maxres),
+ & num_cont(maxres),jcont(maxconts,maxres)
+C 12/26/95 - H-bonding contacts
+ common /contacts_hb/
+ & gacontp_hb1(3,maxconts,maxres),gacontp_hb2(3,maxconts,maxres),
+ & gacontp_hb3(3,maxconts,maxres),
+ & gacontm_hb1(3,maxconts,maxres),gacontm_hb2(3,maxconts,maxres),
+ & gacontm_hb3(3,maxconts,maxres),
+ & gacont_hbr(3,maxconts,maxres),
+ & grij_hb_cont(3,maxconts,maxres),
+ & facont_hb(maxconts,maxres),ees0p(maxconts,maxres),
+ & ees0m(maxconts,maxres),d_cont(maxconts,maxres),
+ & num_cont_hb(maxres),jcont_hb(maxconts,maxres)
+C 9/23/99 Added improper rotation matrices and matrices of dipole-dipole
+C interactions
+C Interactions of pseudo-dipoles generated by loc-el interactions.
+ double precision dip,dipderg,dipderx
+ common /dipint/ dip(4,maxconts,maxres),dipderg(4,maxconts,maxres),
+ & dipderx(3,5,4,maxconts,maxres)
+C 10/30/99 Added other pre-computed vectors and matrices needed
+C to calculate three - six-order el-loc correlation terms
+ double precision Ug,Ugder,Ug2,Ug2der,obrot,obrot2,obrot_der,
+ & obrot2_der,Ub2,Ub2der,mu,muder,EUg,EUgder,CUg,CUgder,
+ & DUg,DUgder,DtUg2,DtUg2der,Ctobr,Ctobrder,Dtobr2,Dtobr2der
+ common /rotat/ Ug(2,2,maxres),Ugder(2,2,maxres),Ug2(2,2,maxres),
+ & Ug2der(2,2,maxres),obrot(2,maxres),obrot2(2,maxres),
+ & obrot_der(2,maxres),obrot2_der(2,maxres)
+C This common block contains vectors and matrices dependent on a single
+C amino-acid residue.
+ common /precomp1/ Ub2(2,maxres),Ub2der(2,maxres),mu(2,maxres),
+ & EUg(2,2,maxres),EUgder(2,2,maxres),CUg(2,2,maxres),
+ & CUgder(2,2,maxres),DUg(2,2,maxres),Dugder(2,2,maxres),
+ & DtUg2(2,2,maxres),DtUg2der(2,2,maxres),Ctobr(2,maxres),
+ & Ctobrder(2,maxres),Dtobr2(2,maxres),Dtobr2der(2,maxres)
+C This common block contains vectors and matrices dependent on two
+C consecutive amino-acid residues.
+ double precision Ug2Db1t,Ug2Db1tder,CUgb2,CUgb2der,EUgC,
+ & EUgCder,EUgD,EUgDder,DtUg2EUg,DtUg2EUgder
+ common /precomp2/ Ug2Db1t(2,maxres),Ug2Db1tder(2,maxres),
+ & CUgb2(2,maxres),CUgb2der(2,maxres),EUgC(2,2,maxres),
+ & EUgCder(2,2,maxres),EUgD(2,2,maxres),EUgDder(2,2,maxres),
+ & DtUg2EUg(2,2,maxres),DtUg2EUgder(2,2,2,maxres),
+ & Ug2DtEUg(2,2,maxres),Ug2DtEUgder(2,2,2,maxres)
+ double precision costab,sintab,costab2,sintab2
+ common /rotat_old/ costab(maxres),sintab(maxres),
+ & costab2(maxres),sintab2(maxres),muder(2,maxres)
+C This common block contains dipole-interaction matrices and their
+C Cartesian derivatives.
+ double precision a_chuj,a_chuj_der
+ common /dipmat/ a_chuj(2,2,maxconts,maxres),
+ & a_chuj_der(2,2,3,5,maxconts,maxres)
+ double precision AEA,AEAderg,AEAderx,AECA,AECAderg,AECAderx,
+ & ADtEA,ADtEAderg,ADtEAderx,AEAb1,AEAb1derg,AEAb1derx,
+ & AEAb2,AEAb2derg,AEAb2derx
+ common /diploc/ AEA(2,2,2),AEAderg(2,2,2),AEAderx(2,2,3,5,2,2),
+ & EAEA(2,2,2), EAEAderg(2,2,2,2), EAEAderx(2,2,3,5,2,2),
+ & AECA(2,2,2),AECAderg(2,2,2),AECAderx(2,2,3,5,2,2),
+ & ADtEA(2,2,2),ADtEAderg(2,2,2,2),ADtEAderx(2,2,3,5,2,2),
+ & ADtEA1(2,2,2),ADtEA1derg(2,2,2,2),ADtEA1derx(2,2,3,5,2,2),
+ & AEAb1(2,2,2),AEAb1derg(2,2,2),AEAb1derx(2,3,5,2,2,2),
+ & AEAb2(2,2,2),AEAb2derg(2,2,2,2),AEAb2derx(2,3,5,2,2,2),
+ & g_contij(3,2),ekont
--- /dev/null
+ double precision sig_comp,chi_comp,chip_comp,sc_cutoff
+ common /contpar/ sig_comp(ntyp,ntyp),chi_comp(ntyp,ntyp),
+ & chip_comp(ntyp,ntyp),sc_cutoff(ntyp,ntyp)
--- /dev/null
+ double precision dcdv,dxdv,dxds,gradx,gradc,gvdwc,gelc,gvdwpp,
+ & gradx_scp,gvdwc_scp,ghpbx,ghpbc,gloc,gvdwx,gradcorr,gradxorr,
+ & gradcorr5,gradcorr6,gel_loc,gcorr3_turn,gcorr4_turn,gcorr6_turn,
+ & gel_loc_loc,gel_loc_turn3,gel_loc_turn4,gel_loc_turn6,gcorr_loc,
+ & g_corr5_loc,g_corr6_loc,gradb,gradbx,gsccorc,gsccorx,gsccor_loc,
+ & gscloc,gsclocx
+ integer nfl,icg
+ logical calc_grad
+ common /derivat/ dcdv(6,maxdim),dxdv(6,maxdim),dxds(6,maxres),
+ & gradx(3,maxres,2),gradc(3,maxres,2),gvdwx(3,maxres),
+ & gvdwc(3,maxres),gelc(3,maxres),gvdwpp(3,maxres),
+ & gradx_scp(3,maxres),
+ & gvdwc_scp(3,maxres),ghpbx(3,maxres),ghpbc(3,maxres),
+ & gloc(maxvar,2),gradcorr(3,maxres),gradxorr(3,maxres),
+ & gradcorr5(3,maxres),gradcorr6(3,maxres),
+ & gel_loc(3,maxres),gcorr3_turn(3,maxres),gcorr4_turn(3,maxres),
+ & gcorr6_turn(3,maxres),gradb(3,maxres),gradbx(3,maxres),
+ & gel_loc_loc(maxvar),gel_loc_turn3(maxvar),gel_loc_turn4(maxvar),
+ & gel_loc_turn6(maxvar),gcorr_loc(maxvar),
+ & g_corr5_loc(maxvar),g_corr6_loc(maxvar),gsccorc(3,maxres),
+ & gsccorx(3,maxres),gsccor_loc(maxres),
+ & gscloc(3,maxres),gsclocx(3,maxres),nfl,icg,calc_grad
+ double precision derx,derx_turn
+ common /deriv_loc/ derx(3,5,2),derx_turn(3,5,2)
+ double precision dXX_C1tab(3,maxres),dYY_C1tab(3,maxres),
+ & dZZ_C1tab(3,maxres),dXX_Ctab(3,maxres),dYY_Ctab(3,maxres),
+ & dZZ_Ctab(3,maxres),dXX_XYZtab(3,maxres),dYY_XYZtab(3,maxres),
+ & dZZ_XYZtab(3,maxres)
+ common /deriv_scloc/ dXX_C1tab,dYY_C1tab,dZZ_C1tab,dXX_Ctab,
+ & dYY_Ctab,dZZ_Ctab,dXX_XYZtab,dYY_XYZtab,dZZ_XYZtab
--- /dev/null
+C-----------------------------------------------------------------------
+C The following COMMON block selects the type of the force field used in
+C calculations and defines weights of various energy terms.
+C 12/1/95 wcorr added
+C-----------------------------------------------------------------------
+ double precision wsc,wscp,welec,wstrain,wtor,wtor_d,wang,wscloc,
+ & wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,wturn3,wturn4,
+ & wturn6,wvdwpp,wbond,weights,scal14,cutoff_corr,delt_corr,
+ & r0_corr
+ integer ipot,n_ene_comp
+ common /ffield/ wsc,wscp,welec,wstrain,wtor,wtor_d,wang,wscloc,
+ & wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,wturn3,wturn4,
+ & wturn6,wvdwpp,wbond,weights(max_ene),
+ & scal14,cutoff_corr,delt_corr,r0_corr,ipot,n_ene_comp
+ common /potentials/ potname(5)
+ character*3 potname
+C-----------------------------------------------------------------------
+C wlong,welec,wtor,wang,wscloc are the weight of the energy terms
+C corresponding to side-chain, electrostatic, torsional, valence-angle,
+C and local side-chain terms.
+C
+C IPOT determines which SC...SC interaction potential will be used:
+C 1 - LJ: 2n-n Lennard-Jones
+C 2 - LJK: 2n-n Kihara type (shifted Lennard-Jones)
+C 3 - BP; Berne-Pechukas (angular dependence)
+C 4 - GB; Gay-Berne (angular dependence)
+C 5 - GBV; Gay-Berne-Vorobjev; angularly-dependent Kihara potential
+C------------------------------------------------------------------------
--- /dev/null
+ integer nbfrag,bfrag,nhfrag,hfrag,bvar_frag,hvar_frag,nhpb0,
+ & nh310frag,h310frag
+ COMMON /c_frag/ nbfrag,bfrag(4,maxres/3),nhfrag,hfrag(2,maxres/3),
+ & nh310frag,h310frag(2,maxres/2)
+ COMMON /frag/ bvar_frag(mxio,6),hvar_frag(mxio,3)
--- /dev/null
+ double precision pi,dwapi,pipol,pi3,dwapi3,deg2rad,rad2deg,angmin
+ common /geo/ pi,dwapi,pipol,pi3,dwapi3,deg2rad,rad2deg,angmin
--- /dev/null
+ character*80 titel
+ common /header/ titel
--- /dev/null
+ double precision aa,bb,augm,aad,bad,app,bpp,ael6,ael3
+ integer nnt,nct,nint_gr,istart,iend,itype,itel,itypro,ielstart,
+ & ielend,nscp_gr,iscpstart,iscpend,iatsc_s,iatsc_e,iatel_s,
+ & iatel_e,iatscp_s,iatscp_e,ispp,iscp,expon,expon2
+ common /interact/aa(ntyp,ntyp),bb(ntyp,ntyp),augm(ntyp,ntyp),
+ & aad(ntyp,2),bad(ntyp,2),app(2,2),bpp(2,2),ael6(2,2),ael3(2,2),
+ & expon,expon2,nnt,nct,nint_gr(maxres),istart(maxres,maxint_gr),
+ & iend(maxres,maxint_gr),itype(maxres),itel(maxres),itypro,
+ & ielstart(maxres),ielend(maxres),nscp_gr(maxres),
+ & iscpstart(maxres,maxint_gr),iscpend(maxres,maxint_gr),
+ & iatsc_s,iatsc_e,iatel_s,iatel_e,iatscp_s,iatscp_e,ispp,iscp
+C 12/1/95 Array EPS included in the COMMON block.
+ double precision eps,sigma,sigmaii,rs0,chi,chip,chip0,alp,signa0,
+ & sigii,sigma0,rr0,r0,r0e,r0d,rpp,epp,elpp6,elpp3,eps_scp,rscp,
+ & eps_orig
+ common /body/eps(ntyp,ntyp),sigma(ntyp,ntyp),sigmaii(ntyp,ntyp),
+ & rs0(ntyp,ntyp),chi(ntyp,ntyp),chip(ntyp),chip0(ntyp),alp(ntyp),
+ & sigma0(ntyp),sigii(ntyp),rr0(ntyp),r0(ntyp,ntyp),r0e(ntyp,ntyp),
+ & r0d(ntyp,2),rpp(2,2),epp(2,2),elpp6(2,2),elpp3(2,2),
+ & eps_scp(20,2),rscp(20,2),eps_orig(ntyp,ntyp)
+c 12/5/03 modified 09/18/03 Bond stretching parameters.
+ double precision vbldp0,vbldsc0,akp,aksc,abond0
+ integer nbondterm
+ common /stretch/ vbldp0,vbldsc0(maxbondterm,ntyp),akp,
+ & aksc(maxbondterm,ntyp),abond0(maxbondterm,ntyp),nbondterm(ntyp)
--- /dev/null
+ double precision a0thet,athet,bthet,polthet,gthet,theta0,sig0,
+ & sigc0,dsc,dsc_inv,bsc,censc,gaussc,dsc0,vbl,vblinv,vblinv2,
+ & vbl_cis,vbl0,vbld_inv
+ integer nlob,loc_start,loc_end,ithet_start,ithet_end,
+ & iphi_start,iphi_end
+C Parameters of the virtual-bond-angle probability distribution
+ common /thetas/ a0thet(ntyp),athet(2,ntyp),bthet(2,ntyp),
+ & polthet(0:3,ntyp),gthet(3,ntyp),theta0(ntyp),sig0(ntyp),
+ & sigc0(ntyp)
+C Parameters of ab initio-derived potential of virtual-bond-angle bending
+ integer nthetyp,ntheterm,ntheterm2,ntheterm3,nsingle,ndouble,
+ & ithetyp(ntyp1),nntheterm
+ double precision aa0thet(maxthetyp1,maxthetyp1,maxthetyp1),
+ & aathet(maxtheterm,maxthetyp1,maxthetyp1,maxthetyp1),
+ & bbthet(maxsingle,maxtheterm2,maxthetyp1,maxthetyp1,maxthetyp1),
+ & ccthet(maxsingle,maxtheterm2,maxthetyp1,maxthetyp1,maxthetyp1),
+ & ddthet(maxsingle,maxtheterm2,maxthetyp1,maxthetyp1,maxthetyp1),
+ & eethet(maxsingle,maxtheterm2,maxthetyp1,maxthetyp1,maxthetyp1),
+ & ffthet(maxdouble,maxdouble,maxtheterm3,maxthetyp1,maxthetyp1,
+ & maxthetyp1),
+ & ggthet(maxdouble,maxdouble,maxtheterm3,maxthetyp1,maxthetyp1,
+ & maxthetyp1)
+ common /theta_abinitio/aa0thet,aathet,bbthet,ccthet,ddthet,eethet,
+ & ffthet,
+ & ggthet,ithetyp,nthetyp,ntheterm,ntheterm2,ntheterm3,nsingle,
+ & ndouble,nntheterm
+C Parameters of the side-chain probability distribution
+ common /sclocal/ dsc(ntyp1),dsc_inv(ntyp1),bsc(maxlob,ntyp),
+ & censc(3,maxlob,ntyp),gaussc(3,3,maxlob,ntyp),dsc0(ntyp1),
+ & nlob(ntyp1)
+C Virtual-bond lenghts
+ common /peptbond/ vbl,vblinv,vblinv2,vbl_cis,vbl0
+ common /indices/ loc_start,loc_end,ithet_start,ithet_end,
+ & iphi_start,iphi_end
+C Inverses of the actual virtual bond lengths
+ common /invlen/ vbld_inv(maxres2)
--- /dev/null
+ double precision tolf,rtolf
+ integer maxfun,maxmin
+ common /minimm/ tolf,rtolf,maxfun,maxmin
--- /dev/null
+ character*3 restyp
+ character*1 onelet
+ common /names/ restyp(ntyp+1),onelet(ntyp+1)
+ character*10 ename,wname
+ integer nprint_ene,print_order
+ common /namterm/ ename(max_ene),wname(max_ene),nprint_ene,
+ & print_order(max_ene)
--- /dev/null
+ double precision ebr,d0cm,akcm,akth,akct,v1ss,v2ss,v3ss,dhpb,
+ & dhpb1,forcon,weidis
+ integer ns,nss,nfree,iss,ihpb,jhpb,nhpb,link_start,link_end,
+ & ibecarb
+ common /sbridge/ ebr,d0cm,akcm,akth,akct,v1ss,v2ss,v3ss,ns,nss,
+ & nfree,iss(maxss)
+ common /links/ dhpb(maxdim),dhpb1(maxdim),forcon(maxdim),
+ & ihpb(maxdim),jhpb(maxdim),ibecarb(maxdim),nhpb
+ common /restraints/ weidis
+ common /links_split/ link_start,link_end
--- /dev/null
+C Parameters of the SCCOR term
+ double precision v1sccor,v2sccor
+ integer nterm_sccor
+ common/torsion/v1sccor(maxterm_sccor,20,20),
+ & v2sccor(maxterm_sccor,20,20),
+ & nterm_sccor
--- /dev/null
+C Parameters of the SC rotamers (local) term
+ double precision sc_parmin
+ common/scrot/sc_parmin(maxsccoef,20)
--- /dev/null
+ DOUBLE PRECISION BATIME,TIMLIM,STIME,PREVTIM,SAFETY,RSTIME
+ INTEGER WhatsUp,ndelta
+ logical cutoffviol,cutoffeval,llocal
+ COMMON/TIME1/STIME,TIMLIM,BATIME,PREVTIM,SAFETY,RSTIME
+ COMMON/STOPTIM/WhatsUp,ndelta,cutoffviol,cutoffeval,llocal
+ double precision t_func,t_grad,t_fhel,t_fbet,t_ghel,t_gbet,t_viol,
+ & t_gviol,t_map,t_alamap,t_betamap
+ integer n_func,n_grad,n_fhel,n_fbet,n_ghel,n_gbet,n_viol,n_gviol,
+ & n_map,n_alamap,n_betamap
+ common /timing/ t_func,t_grad,t_fhel,t_fbet,t_ghel,t_gbet,t_viol,
+ & t_gviol,t_map,t_alamap,t_betamap,
+ & n_func,n_grad,n_fhel,n_fbet,n_ghel,n_gbet,n_viol,n_gviol,
+ & n_map,n_alamap,n_betamap
--- /dev/null
+ integer ndih_constr,idih_constr(maxdih_constr)
+ integer ndih_nconstr,idih_nconstr(maxdih_constr)
+ double precision phi0(maxdih_constr),drange(maxdih_constr),ftors
+ common /torcnstr/ phi0,drange,ftors,ndih_constr,idih_constr,
+ & ndih_nconstr,idih_nconstr
--- /dev/null
+C Torsional constants of the rotation about virtual-bond dihedral angles
+ double precision v1,v2,vlor1,vlor2,vlor3,v0
+ integer itortyp,ntortyp,nterm,nlor,nterm_old
+ common/torsion/v0(maxtor,maxtor),v1(maxterm,maxtor,maxtor),
+ & v2(maxterm,maxtor,maxtor),vlor1(maxlor,maxtor,maxtor),
+ & vlor2(maxlor,maxtor,maxtor),vlor3(maxlor,maxtor,maxtor),
+ & itortyp(ntyp),ntortyp,nterm(maxtor,maxtor),nlor(maxtor,maxtor)
+ & ,nterm_old
+C 6/23/01 - constants for double torsionals
+ double precision v1c,v1s,v2c,v2s
+ integer ntermd_1,ntermd_2
+ common /torsiond/ v1c(2,maxtermd_1,maxtor,maxtor,maxtor),
+ & v1s(2,maxtermd_1,maxtor,maxtor,maxtor),
+ & v2c(maxtermd_2,maxtermd_2,maxtor,maxtor,maxtor),
+ & v2s(maxtermd_2,maxtermd_2,maxtor,maxtor,maxtor),
+ & ntermd_1(maxtor,maxtor,maxtor),ntermd_2(maxtor,maxtor,maxtor)
+C 9/18/99 - added Fourier coeffficients of the expansion of local energy
+C surface
+ double precision b1,b2,cc,dd,ee,ctilde,dtilde,b1tilde
+ integer nloctyp
+ common/fourier/ b1(2,maxtor),b2(2,maxtor),cc(2,2,maxtor),
+ & dd(2,2,maxtor),ee(2,2,maxtor),ctilde(2,2,maxtor),
+ & dtilde(2,2,maxtor),b1tilde(2,maxtor),nloctyp
+ double precision b
+ common /fourier1/ b(13,maxtor)
--- /dev/null
+ common /vectors/ uy(3,maxres),uz(3,maxres),
+ & uygrad(3,3,2,maxres),uzgrad(3,3,2,maxres)
+
--- /dev/null
+ double precision ww,ww0,ww_low,ww_up,ww_orig,x_orig,
+ & epp_low,epp_up,rpp_low,rpp_up,elpp6_low,elpp6_up,elpp3_low,
+ & elpp3_up,b_low,b_up,epscp_low,epscp_up,rscp_low,rscp_up,
+ & x_up,x_low,xm,xm1,xm2,epss_low,epss_up,epsp_low,epsp_up
+ integer imask,mask_elec,mask_fourier,mod_fourier,mask_scp,indz,iw,
+ & nsingle_sc,npair_sc,ityp_ssc,ityp_psc
+ logical mod_other_params,mod_elec,mod_scp,mod_side
+ common /chujec/ ww(max_ene),ww0(max_ene),ww_low(max_ene),
+ & ww_up(max_ene),ww_orig(max_ene),x_orig(max_paropt),
+ & epp_low(2,2),epp_up(2,2),rpp_low(2,2),rpp_up(2,2),
+ & elpp6_low(2,2),elpp6_up(2,2),elpp3_low(2,2),elpp3_up(2,2),
+ & b_low(13,3),b_up(13,3),x_up(max_paropt),x_low(max_paropt),
+ & epscp_low(0:20,2),epscp_up(0:20,2),rscp_low(0:20,2),
+ & rscp_up(0:20,2),epss_low(ntyp),epss_up(ntyp),epsp_low(nntyp),
+ & epsp_up(nntyp),
+ & xm(max_paropt,0:maxprot),xm1(max_paropt,0:maxprot),
+ & xm2(max_paropt,0:maxprot),
+ & imask(max_ene),nsingle_sc,npair_sc,ityp_ssc(ntyp),
+ & ityp_psc(2,nntyp),mask_elec(2,2,4),
+ & mask_fourier(13,3),
+ & mask_scp(0:20,2,2),mod_other_params,mod_fourier(0:3),
+ & mod_elec,mod_scp,mod_side,indz(maxbatch+1,maxprot),iw(max_ene)
call etotal(energia(0),fT)
totfree(i)=energia(0)
#ifdef DEBUG
- write (iout,*) i," energia",(energia(j),j=0,19)
+ write (iout,*) i," energia",(energia(j),j=0,21)
+ call enerprint(energia(0),ft)
+ call flush(iout)
#endif
do k=1,max_ene
enetb(k,i)=energia(k)
write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
lprint_cart=index(controlcard,"PRINT_CART") .gt.0
lprint_int=index(controlcard,"PRINT_INT") .gt.0
+ with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
+ call readi(controlcard,'CONSTR_DIST',constr_dist,0)
+ write (iout,*) "with_dihed_constr ",with_dihed_constr,
+ & " CONSTR_DIST",constr_dist
+ call flush(iout)
if (min_var) iopt=1
return
end
include 'COMMON.CONTROL'
include 'COMMON.CONTACTS'
include 'COMMON.TIME1'
+ include 'COMMON.TORCNSTR'
#ifdef MPL
include 'COMMON.INFO'
#endif
print *,'Call Read_Bridge.'
call read_bridge
+ if (with_dihed_constr) then
+
+ read (inp,*) ndih_constr
+ if (ndih_constr.gt.0) then
+ read (inp,*) ftors
+ write (iout,*) 'FTORS',ftors
+ read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
+ write (iout,*)
+ & 'There are',ndih_constr,' constraints on phi angles.'
+ do i=1,ndih_constr
+ write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
+ enddo
+ do i=1,ndih_constr
+ phi0(i)=deg2rad*phi0(i)
+ drange(i)=deg2rad*drange(i)
+ enddo
+ endif
+
+ endif
nnt=1
nct=nres
print *,'NNT=',NNT,' NCT=',NCT
endif
call contact(.true.,ncont_ref,icont_ref)
endif
+c Read distance restraints
+ if (constr_dist.gt.0) then
+ call read_dist_constr
+ call hpb_partition
+ endif
return
end
c-----------------------------------------------------------------------------
enddo
write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
20 continue
- dhpb(i)=dbr
- forcon(i)=fbr
+c dhpb(i)=dbr
+c forcon(i)=fbr
enddo
do i=1,nss
ihpb(i)=ihpb(i)+nres
return
end
c----------------------------------------------------------------------------
+ subroutine multreadi(rekord,lancuch,tablica,dim,default)
+ implicit none
+ integer dim,i
+ integer tablica(dim),default
+ character*(*) rekord,lancuch
+ character*80 aux
+ integer ilen,iread
+ external ilen
+ do i=1,dim
+ tablica(i)=default
+ enddo
+ iread=index(rekord,lancuch(:ilen(lancuch))//"=")
+ if (iread.eq.0) return
+ iread=iread+ilen(lancuch)+1
+ read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
+ 10 return
+ end
+c----------------------------------------------------------------------------
subroutine card_concat(card)
include 'DIMENSIONS'
include 'COMMON.IOUNITS'
#endif
return
end
+c-------------------------------------------------------------------------------
+ subroutine read_dist_constr
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CONTROL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.SBRIDGE'
+ integer ifrag_(2,100),ipair_(2,100)
+ double precision wfrag_(100),wpair_(100)
+ character*500 controlcard
+c write (iout,*) "Calling read_dist_constr"
+c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
+ call card_concat(controlcard)
+c call flush(iout)
+c write (iout,'(a)') controlcard
+ call readi(controlcard,"NFRAG",nfrag_,0)
+ call readi(controlcard,"NPAIR",npair_,0)
+ call readi(controlcard,"NDIST",ndist_,0)
+ call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
+ call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
+ call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
+ call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
+ call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
+ write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
+ write (iout,*) "IFRAG"
+ do i=1,nfrag_
+ write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
+ enddo
+ write (iout,*) "IPAIR"
+ do i=1,npair_
+ write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
+ enddo
+ call flush(iout)
+ if (.not.refstr .and. nfrag_.gt.0) then
+ write (iout,*)
+ & "ERROR: no reference structure to compute distance restraints"
+ write (iout,*)
+ & "Restraints must be specified explicitly (NDIST=number)"
+ stop
+ endif
+ if (nfrag_.lt.2 .and. npair_.gt.0) then
+ write (iout,*) "ERROR: Less than 2 fragments specified",
+ & " but distance restraints between pairs requested"
+ stop
+ endif
+ call flush(iout)
+ do i=1,nfrag_
+ if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
+ if (ifrag_(2,i).gt.nstart_sup+nsup-1)
+ & ifrag_(2,i)=nstart_sup+nsup-1
+c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
+ call flush(iout)
+ if (wfrag_(i).gt.0.0d0) then
+ do j=ifrag_(1,i),ifrag_(2,i)-1
+ do k=j+1,ifrag_(2,i)
+ write (iout,*) "j",j," k",k
+ ddjk=dist(j,k)
+ if (constr_dist.eq.1) then
+ nhpb=nhpb+1
+ ihpb(nhpb)=j
+ jhpb(nhpb)=k
+ dhpb(nhpb)=ddjk
+ forcon(nhpb)=wfrag_(i)
+ else if (constr_dist.eq.2) then
+ if (ddjk.le.dist_cut) then
+ nhpb=nhpb+1
+ ihpb(nhpb)=j
+ jhpb(nhpb)=k
+ dhpb(nhpb)=ddjk
+ forcon(nhpb)=wfrag_(i)
+ endif
+ else
+ nhpb=nhpb+1
+ ihpb(nhpb)=j
+ jhpb(nhpb)=k
+ dhpb(nhpb)=ddjk
+ forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
+ endif
+ write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
+ & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+ enddo
+ enddo
+ endif
+ enddo
+ do i=1,npair_
+ if (wpair_(i).gt.0.0d0) then
+ ii = ipair_(1,i)
+ jj = ipair_(2,i)
+ if (ii.gt.jj) then
+ itemp=ii
+ ii=jj
+ jj=itemp
+ endif
+ do j=ifrag_(1,ii),ifrag_(2,ii)
+ do k=ifrag_(1,jj),ifrag_(2,jj)
+ nhpb=nhpb+1
+ ihpb(nhpb)=j
+ jhpb(nhpb)=k
+ forcon(nhpb)=wpair_(i)
+ dhpb(nhpb)=dist(j,k)
+ write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
+ & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+ enddo
+ enddo
+ endif
+ enddo
+ do i=1,ndist_
+ read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
+ & ibecarb(i),forcon(nhpb+1)
+ if (forcon(nhpb+1).gt.0.0d0) then
+ nhpb=nhpb+1
+ if (ibecarb(i).gt.0) then
+ ihpb(i)=ihpb(i)+nres
+ jhpb(i)=jhpb(i)+nres
+ endif
+ if (dhpb(nhpb).eq.0.0d0)
+ & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
+ endif
+ enddo
+ do i=1,nhpb
+ write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
+ & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
+ enddo
+ call flush(iout)
+ return
+ end
--- /dev/null
+# This make file is part of the xdrf package.
+#
+# (C) 1995 Frans van Hoesel, hoesel@chem.rug.nl
+#
+# 2006 modified by Cezary Czaplewski
+
+# Set C compiler and flags for ARCH
+CC = gcc
+CFLAGS = -O
+
+M4 = m4
+M4FILE = underscore.m4
+
+libxdrf.a: libxdrf.o ftocstr.o
+ ar cr libxdrf.a $?
+
+clean:
+ rm -f libxdrf.o ftocstr.o libxdrf.a
+
+ftocstr.o: ftocstr.c
+ $(CC) $(CFLAGS) -c ftocstr.c
+
+libxdrf.o: libxdrf.m4 $(M4FILE)
+ $(M4) $(M4FILE) libxdrf.m4 > libxdrf.c
+ $(CC) $(CFLAGS) -c libxdrf.c
+ rm -f libxdrf.c
+
--- /dev/null
+# This make file is part of the xdrf package.
+#
+# (C) 1995 Frans van Hoesel, hoesel@chem.rug.nl
+#
+# 2006 modified by Cezary Czaplewski
+
+# Set C compiler and flags for ARCH
+BGLSYS = /bgl/BlueLight/ppcfloor/bglsys
+
+CC = /usr/bin/blrts_xlc
+CPPC = /usr/bin/blrts_xlc
+
+CFLAGS= -O2 -I$(BGLSYS)/include -L$(BGLSYS)/lib -qarch=440d -qtune=440
+
+M4 = m4
+M4FILE = RS6K.m4
+
+libxdrf.a: libxdrf.o ftocstr.o xdr_array.o xdr.o xdr_float.o xdr_stdio.o
+ ar cr libxdrf.a $?
+
+clean:
+ rm -f *.o libxdrf.a
+
+ftocstr.o: ftocstr.c
+ $(CC) $(CFLAGS) -c ftocstr.c
+
+libxdrf.o: libxdrf.m4 $(M4FILE)
+ $(M4) $(M4FILE) libxdrf.m4 > libxdrf.c
+ $(CC) $(CFLAGS) -c libxdrf.c
+# rm -f libxdrf.c
+
--- /dev/null
+# This make file is part of the xdrf package.
+#
+# (C) 1995 Frans van Hoesel, hoesel@chem.rug.nl
+#
+# 2006 modified by Cezary Czaplewski
+
+# Set C compiler and flags for ARCH
+CC = cc
+CFLAGS = -O
+
+M4 = m4
+M4FILE = underscore.m4
+
+libxdrf.a: libxdrf.o ftocstr.o
+ ar cr libxdrf.a $?
+
+clean:
+ rm -f libxdrf.o ftocstr.o libxdrf.a
+
+ftocstr.o: ftocstr.c
+ $(CC) $(CFLAGS) -c ftocstr.c
+
+libxdrf.o: libxdrf.m4 $(M4FILE)
+ $(M4) $(M4FILE) libxdrf.m4 > libxdrf.c
+ $(CC) $(CFLAGS) -c libxdrf.c
+ rm -f libxdrf.c
+
--- /dev/null
+divert(-1)
+undefine(`len')
+#
+# do nothing special to FORTRAN function names
+#
+define(`FUNCTION',`$1')
+#
+# FORTRAN character strings are passed as follows:
+# a pointer to the base of the string is passed in the normal
+# argument list, and the length is passed by value as an extra
+# argument, after all of the other arguments.
+#
+define(`ARGS',`($1`'undivert(1))')
+define(`SAVE',`divert(1)$1`'divert(0)')
+define(`STRING_ARG',`$1_ptr`'SAVE(`, $1_len')')
+define(`STRING_ARG_DECL',`char * $1_ptr; int $1_len')
+define(`STRING_LEN',`$1_len')
+define(`STRING_PTR',`$1_ptr')
+divert(0)
+
--- /dev/null
+
+
+int ftocstr(ds, dl, ss, sl)
+ char *ds, *ss; /* dst, src ptrs */
+ int dl; /* dst max len */
+ int sl; /* src len */
+{
+ char *p;
+
+ for (p = ss + sl; --p >= ss && *p == ' '; ) ;
+ sl = p - ss + 1;
+ dl--;
+ ds[0] = 0;
+ if (sl > dl)
+ return 1;
+ while (sl--)
+ (*ds++ = *ss++);
+ *ds = '\0';
+ return 0;
+}
+
+
+int ctofstr(ds, dl, ss)
+ char *ds; /* dest space */
+ int dl; /* max dest length */
+ char *ss; /* src string (0-term) */
+{
+ while (dl && *ss) {
+ *ds++ = *ss++;
+ dl--;
+ }
+ while (dl--)
+ *ds++ = ' ';
+ return 0;
+}
--- /dev/null
+/*____________________________________________________________________________
+ |
+ | libxdrf - portable fortran interface to xdr. some xdr routines
+ | are C routines for compressed coordinates
+ |
+ | version 1.1
+ |
+ | This collection of routines is intended to write and read
+ | data in a portable way to a file, so data written on one type
+ | of machine can be read back on a different type.
+ |
+ | all fortran routines use an integer 'xdrid', which is an id to the
+ | current xdr file, and is set by xdrfopen.
+ | most routines have in integer 'ret' which is the return value.
+ | The value of 'ret' is zero on failure, and most of the time one
+ | on succes.
+ |
+ | There are three routines useful for C users:
+ | xdropen(), xdrclose(), xdr3dfcoord().
+ | The first two replace xdrstdio_create and xdr_destroy, and *must* be
+ | used when you plan to use xdr3dfcoord(). (they are also a bit
+ | easier to interface). For writing data other than compressed coordinates
+ | you should use the standard C xdr routines (see xdr man page)
+ |
+ | xdrfopen(xdrid, filename, mode, ret)
+ | character *(*) filename
+ | character *(*) mode
+ |
+ | this will open the file with the given filename (string)
+ | and the given mode, it returns an id in xdrid, which is
+ | to be used in all other calls to xdrf routines.
+ | mode is 'w' to create, or update an file, for all other
+ | values of mode the file is opened for reading
+ |
+ | you need to call xdrfclose to flush the output and close
+ | the file.
+ | Note that you should not use xdrstdio_create, which comes with the
+ | standard xdr library
+ |
+ | xdrfclose(xdrid, ret)
+ | flush the data to the file, and closes the file;
+ | You should not use xdr_destroy (which comes standard with
+ | the xdr libraries.
+ |
+ | xdrfbool(xdrid, bp, ret)
+ | integer pb
+ |
+ | This filter produces values of either 1 or 0
+ |
+ | xdrfchar(xdrid, cp, ret)
+ | character cp
+ |
+ | filter that translate between characters and their xdr representation
+ | Note that the characters in not compressed and occupies 4 bytes.
+ |
+ | xdrfdouble(xdrid, dp, ret)
+ | double dp
+ |
+ | read/write a double.
+ |
+ | xdrffloat(xdrid, fp, ret)
+ | float fp
+ |
+ | read/write a float.
+ |
+ | xdrfint(xdrid, ip, ret)
+ | integer ip
+ |
+ | read/write integer.
+ |
+ | xdrflong(xdrid, lp, ret)
+ | integer lp
+ |
+ | this routine has a possible portablility problem due to 64 bits longs.
+ |
+ | xdrfshort(xdrid, sp, ret)
+ | integer *2 sp
+ |
+ | xdrfstring(xdrid, sp, maxsize, ret)
+ | character *(*)
+ | integer maxsize
+ |
+ | read/write a string, with maximum length given by maxsize
+ |
+ | xdrfwrapstring(xdris, sp, ret)
+ | character *(*)
+ |
+ | read/write a string (it is the same as xdrfstring accept that it finds
+ | the stringlength itself.
+ |
+ | xdrfvector(xdrid, cp, size, xdrfproc, ret)
+ | character *(*)
+ | integer size
+ | external xdrfproc
+ |
+ | read/write an array pointed to by cp, with number of elements
+ | defined by 'size'. the routine 'xdrfproc' is the name
+ | of one of the above routines to read/write data (like xdrfdouble)
+ | In contrast with the c-version you don't need to specify the
+ | byte size of an element.
+ | xdrfstring is not allowed here (it is in the c version)
+ |
+ | xdrf3dfcoord(xdrid, fp, size, precision, ret)
+ | real (*) fp
+ | real precision
+ | integer size
+ |
+ | this is *NOT* a standard xdr routine. I named it this way, because
+ | it invites people to use the other xdr routines.
+ | It is introduced to store specifically 3d coordinates of molecules
+ | (as found in molecular dynamics) and it writes it in a compressed way.
+ | It starts by multiplying all numbers by precision and
+ | rounding the result to integer. effectively converting
+ | all floating point numbers to fixed point.
+ | it uses an algorithm for compression that is optimized for
+ | molecular data, but could be used for other 3d coordinates
+ | as well. There is subtantial overhead involved, so call this
+ | routine only if you have a large number of coordinates to read/write
+ |
+ | ________________________________________________________________________
+ |
+ | Below are the routines to be used by C programmers. Use the 'normal'
+ | xdr routines to write integers, floats, etc (see man xdr)
+ |
+ | int xdropen(XDR *xdrs, const char *filename, const char *type)
+ | This will open the file with the given filename and the
+ | given mode. You should pass it an allocated XDR struct
+ | in xdrs, to be used in all other calls to xdr routines.
+ | Mode is 'w' to create, or update an file, and for all
+ | other values of mode the file is opened for reading.
+ | You need to call xdrclose to flush the output and close
+ | the file.
+ |
+ | Note that you should not use xdrstdio_create, which
+ | comes with the standard xdr library.
+ |
+ | int xdrclose(XDR *xdrs)
+ | Flush the data to the file, and close the file;
+ | You should not use xdr_destroy (which comes standard
+ | with the xdr libraries).
+ |
+ | int xdr3dfcoord(XDR *xdrs, float *fp, int *size, float *precision)
+ | This is \fInot\fR a standard xdr routine. I named it this
+ | way, because it invites people to use the other xdr
+ | routines.
+ |
+ | (c) 1995 Frans van Hoesel, hoesel@chem.rug.nl
+*/
+
+
+#include <limits.h>
+#include <malloc.h>
+#include <math.h>
+/* #include <rpc/rpc.h>
+#include <rpc/xdr.h> */
+#include "xdr.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include "xdrf.h"
+
+int ftocstr(char *, int, char *, int);
+int ctofstr(char *, int, char *);
+
+#define MAXID 20
+static FILE *xdrfiles[MAXID];
+static XDR *xdridptr[MAXID];
+static char xdrmodes[MAXID];
+static unsigned int cnt;
+
+typedef void (* FUNCTION(xdrfproc)) (int *, void *, int *);
+
+void
+FUNCTION(xdrfbool) ARGS(`xdrid, pb, ret')
+int *xdrid, *ret;
+int *pb;
+{
+ *ret = xdr_bool(xdridptr[*xdrid], (bool_t *) pb);
+ cnt += sizeof(int);
+}
+
+void
+FUNCTION(xdrfchar) ARGS(`xdrid, cp, ret')
+int *xdrid, *ret;
+char *cp;
+{
+ *ret = xdr_char(xdridptr[*xdrid], cp);
+ cnt += sizeof(char);
+}
+
+void
+FUNCTION(xdrfdouble) ARGS(`xdrid, dp, ret')
+int *xdrid, *ret;
+double *dp;
+{
+ *ret = xdr_double(xdridptr[*xdrid], dp);
+ cnt += sizeof(double);
+}
+
+void
+FUNCTION(xdrffloat) ARGS(`xdrid, fp, ret')
+int *xdrid, *ret;
+float *fp;
+{
+ *ret = xdr_float(xdridptr[*xdrid], fp);
+ cnt += sizeof(float);
+}
+
+void
+FUNCTION(xdrfint) ARGS(`xdrid, ip, ret')
+int *xdrid, *ret;
+int *ip;
+{
+ *ret = xdr_int(xdridptr[*xdrid], ip);
+ cnt += sizeof(int);
+}
+
+void
+FUNCTION(xdrflong) ARGS(`xdrid, lp, ret')
+int *xdrid, *ret;
+long *lp;
+{
+ *ret = xdr_long(xdridptr[*xdrid], lp);
+ cnt += sizeof(long);
+}
+
+void
+FUNCTION(xdrfshort) ARGS(`xdrid, sp, ret')
+int *xdrid, *ret;
+short *sp;
+{
+ *ret = xdr_short(xdridptr[*xdrid], sp);
+ cnt += sizeof(sp);
+}
+
+void
+FUNCTION(xdrfuchar) ARGS(`xdrid, ucp, ret')
+int *xdrid, *ret;
+char *ucp;
+{
+ *ret = xdr_u_char(xdridptr[*xdrid], ucp);
+ cnt += sizeof(char);
+}
+
+void
+FUNCTION(xdrfulong) ARGS(`xdrid, ulp, ret')
+int *xdrid, *ret;
+unsigned long *ulp;
+{
+ *ret = xdr_u_long(xdridptr[*xdrid], ulp);
+ cnt += sizeof(unsigned long);
+}
+
+void
+FUNCTION(xdrfushort) ARGS(`xdrid, usp, ret')
+int *xdrid, *ret;
+unsigned short *usp;
+{
+ *ret = xdr_u_short(xdridptr[*xdrid], usp);
+ cnt += sizeof(unsigned short);
+}
+
+void
+FUNCTION(xdrf3dfcoord) ARGS(`xdrid, fp, size, precision, ret')
+int *xdrid, *ret;
+float *fp;
+int *size;
+float *precision;
+{
+ *ret = xdr3dfcoord(xdridptr[*xdrid], fp, size, precision);
+}
+
+void
+FUNCTION(xdrfstring) ARGS(`xdrid, STRING_ARG(sp), maxsize, ret')
+int *xdrid, *ret;
+STRING_ARG_DECL(sp);
+int *maxsize;
+{
+ char *tsp;
+
+ tsp = (char*) malloc(((STRING_LEN(sp)) + 1) * sizeof(char));
+ if (tsp == NULL) {
+ *ret = -1;
+ return;
+ }
+ if (ftocstr(tsp, *maxsize+1, STRING_PTR(sp), STRING_LEN(sp))) {
+ *ret = -1;
+ free(tsp);
+ return;
+ }
+ *ret = xdr_string(xdridptr[*xdrid], (char **) &tsp, (u_int) *maxsize);
+ ctofstr( STRING_PTR(sp), STRING_LEN(sp), tsp);
+ cnt += *maxsize;
+ free(tsp);
+}
+
+void
+FUNCTION(xdrfwrapstring) ARGS(`xdrid, STRING_ARG(sp), ret')
+int *xdrid, *ret;
+STRING_ARG_DECL(sp);
+{
+ char *tsp;
+ int maxsize;
+ maxsize = (STRING_LEN(sp)) + 1;
+ tsp = (char*) malloc(maxsize * sizeof(char));
+ if (tsp == NULL) {
+ *ret = -1;
+ return;
+ }
+ if (ftocstr(tsp, maxsize, STRING_PTR(sp), STRING_LEN(sp))) {
+ *ret = -1;
+ free(tsp);
+ return;
+ }
+ *ret = xdr_string(xdridptr[*xdrid], (char **) &tsp, (u_int)maxsize);
+ ctofstr( STRING_PTR(sp), STRING_LEN(sp), tsp);
+ cnt += maxsize;
+ free(tsp);
+}
+
+void
+FUNCTION(xdrfopaque) ARGS(`xdrid, cp, ccnt, ret')
+int *xdrid, *ret;
+caddr_t *cp;
+int *ccnt;
+{
+ *ret = xdr_opaque(xdridptr[*xdrid], (caddr_t)*cp, (u_int)*ccnt);
+ cnt += *ccnt;
+}
+
+void
+FUNCTION(xdrfsetpos) ARGS(`xdrid, pos, ret')
+int *xdrid, *ret;
+int *pos;
+{
+ *ret = xdr_setpos(xdridptr[*xdrid], (u_int) *pos);
+}
+
+void
+FUNCTION(xdrf) ARGS(`xdrid, pos')
+int *xdrid, *pos;
+{
+ *pos = xdr_getpos(xdridptr[*xdrid]);
+}
+
+void
+FUNCTION(xdrfvector) ARGS(`xdrid, cp, size, elproc, ret')
+int *xdrid, *ret;
+char *cp;
+int *size;
+FUNCTION(xdrfproc) elproc;
+{
+ int lcnt;
+ cnt = 0;
+ for (lcnt = 0; lcnt < *size; lcnt++) {
+ elproc(xdrid, (cp+cnt) , ret);
+ }
+}
+
+
+void
+FUNCTION(xdrfclose) ARGS(`xdrid, ret')
+int *xdrid;
+int *ret;
+{
+ *ret = xdrclose(xdridptr[*xdrid]);
+ cnt = 0;
+}
+
+void
+FUNCTION(xdrfopen) ARGS(`xdrid, STRING_ARG(fp), STRING_ARG(mode), ret')
+int *xdrid;
+STRING_ARG_DECL(fp);
+STRING_ARG_DECL(mode);
+int *ret;
+{
+ char fname[512];
+ char fmode[3];
+
+ if (ftocstr(fname, sizeof(fname), STRING_PTR(fp), STRING_LEN(fp))) {
+ *ret = 0;
+ }
+ if (ftocstr(fmode, sizeof(fmode), STRING_PTR(mode),
+ STRING_LEN(mode))) {
+ *ret = 0;
+ }
+
+ *xdrid = xdropen(NULL, fname, fmode);
+ if (*xdrid == 0)
+ *ret = 0;
+ else
+ *ret = 1;
+}
+
+/*___________________________________________________________________________
+ |
+ | what follows are the C routines for opening, closing xdr streams
+ | and the routine to read/write compressed coordinates together
+ | with some routines to assist in this task (those are marked
+ | static and cannot be called from user programs)
+*/
+#define MAXABS INT_MAX-2
+
+#ifndef MIN
+#define MIN(x,y) ((x) < (y) ? (x):(y))
+#endif
+#ifndef MAX
+#define MAX(x,y) ((x) > (y) ? (x):(y))
+#endif
+#ifndef SQR
+#define SQR(x) ((x)*(x))
+#endif
+static int magicints[] = {
+ 0, 0, 0, 0, 0, 0, 0, 0, 0,
+ 8, 10, 12, 16, 20, 25, 32, 40, 50, 64,
+ 80, 101, 128, 161, 203, 256, 322, 406, 512, 645,
+ 812, 1024, 1290, 1625, 2048, 2580, 3250, 4096, 5060, 6501,
+ 8192, 10321, 13003, 16384, 20642, 26007, 32768, 41285, 52015, 65536,
+ 82570, 104031, 131072, 165140, 208063, 262144, 330280, 416127, 524287, 660561,
+ 832255, 1048576, 1321122, 1664510, 2097152, 2642245, 3329021, 4194304, 5284491, 6658042,
+ 8388607, 10568983, 13316085, 16777216 };
+
+#define FIRSTIDX 9
+/* note that magicints[FIRSTIDX-1] == 0 */
+#define LASTIDX (sizeof(magicints) / sizeof(*magicints))
+
+
+/*__________________________________________________________________________
+ |
+ | xdropen - open xdr file
+ |
+ | This versions differs from xdrstdio_create, because I need to know
+ | the state of the file (read or write) so I can use xdr3dfcoord
+ | in eigther read or write mode, and the file descriptor
+ | so I can close the file (something xdr_destroy doesn't do).
+ |
+*/
+
+int xdropen(XDR *xdrs, const char *filename, const char *type) {
+ static int init_done = 0;
+ enum xdr_op lmode;
+ const char *type1;
+ int xdrid;
+
+ if (init_done == 0) {
+ for (xdrid = 1; xdrid < MAXID; xdrid++) {
+ xdridptr[xdrid] = NULL;
+ }
+ init_done = 1;
+ }
+ xdrid = 1;
+ while (xdrid < MAXID && xdridptr[xdrid] != NULL) {
+ xdrid++;
+ }
+ if (xdrid == MAXID) {
+ return 0;
+ }
+ if (*type == 'w' || *type == 'W') {
+ type = "w+";
+ type1 = "w+";
+ lmode = XDR_ENCODE;
+ } else if (*type == 'a' || *type == 'A') {
+ type = "w+";
+ type1 = "a+";
+ lmode = XDR_ENCODE;
+ } else {
+ type = "r";
+ type1 = "r";
+ lmode = XDR_DECODE;
+ }
+ xdrfiles[xdrid] = fopen(filename, type1);
+ if (xdrfiles[xdrid] == NULL) {
+ xdrs = NULL;
+ return 0;
+ }
+ xdrmodes[xdrid] = *type;
+ /* next test isn't usefull in the case of C language
+ * but is used for the Fortran interface
+ * (C users are expected to pass the address of an already allocated
+ * XDR staructure)
+ */
+ if (xdrs == NULL) {
+ xdridptr[xdrid] = (XDR *) malloc(sizeof(XDR));
+ xdrstdio_create(xdridptr[xdrid], xdrfiles[xdrid], lmode);
+ } else {
+ xdridptr[xdrid] = xdrs;
+ xdrstdio_create(xdrs, xdrfiles[xdrid], lmode);
+ }
+ return xdrid;
+}
+
+/*_________________________________________________________________________
+ |
+ | xdrclose - close a xdr file
+ |
+ | This will flush the xdr buffers, and destroy the xdr stream.
+ | It also closes the associated file descriptor (this is *not*
+ | done by xdr_destroy).
+ |
+*/
+
+int xdrclose(XDR *xdrs) {
+ int xdrid;
+
+ if (xdrs == NULL) {
+ fprintf(stderr, "xdrclose: passed a NULL pointer\n");
+ exit(1);
+ }
+ for (xdrid = 1; xdrid < MAXID; xdrid++) {
+ if (xdridptr[xdrid] == xdrs) {
+
+ xdr_destroy(xdrs);
+ fclose(xdrfiles[xdrid]);
+ xdridptr[xdrid] = NULL;
+ return 1;
+ }
+ }
+ fprintf(stderr, "xdrclose: no such open xdr file\n");
+ exit(1);
+
+}
+
+/*____________________________________________________________________________
+ |
+ | sendbits - encode num into buf using the specified number of bits
+ |
+ | This routines appends the value of num to the bits already present in
+ | the array buf. You need to give it the number of bits to use and you
+ | better make sure that this number of bits is enough to hold the value
+ | Also num must be positive.
+ |
+*/
+
+static void sendbits(int buf[], int num_of_bits, int num) {
+
+ unsigned int cnt, lastbyte;
+ int lastbits;
+ unsigned char * cbuf;
+
+ cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
+ cnt = (unsigned int) buf[0];
+ lastbits = buf[1];
+ lastbyte =(unsigned int) buf[2];
+ while (num_of_bits >= 8) {
+ lastbyte = (lastbyte << 8) | ((num >> (num_of_bits -8)) /* & 0xff*/);
+ cbuf[cnt++] = lastbyte >> lastbits;
+ num_of_bits -= 8;
+ }
+ if (num_of_bits > 0) {
+ lastbyte = (lastbyte << num_of_bits) | num;
+ lastbits += num_of_bits;
+ if (lastbits >= 8) {
+ lastbits -= 8;
+ cbuf[cnt++] = lastbyte >> lastbits;
+ }
+ }
+ buf[0] = cnt;
+ buf[1] = lastbits;
+ buf[2] = lastbyte;
+ if (lastbits>0) {
+ cbuf[cnt] = lastbyte << (8 - lastbits);
+ }
+}
+
+/*_________________________________________________________________________
+ |
+ | sizeofint - calculate bitsize of an integer
+ |
+ | return the number of bits needed to store an integer with given max size
+ |
+*/
+
+static int sizeofint(const int size) {
+ unsigned int num = 1;
+ int num_of_bits = 0;
+
+ while (size >= num && num_of_bits < 32) {
+ num_of_bits++;
+ num <<= 1;
+ }
+ return num_of_bits;
+}
+
+/*___________________________________________________________________________
+ |
+ | sizeofints - calculate 'bitsize' of compressed ints
+ |
+ | given the number of small unsigned integers and the maximum value
+ | return the number of bits needed to read or write them with the
+ | routines receiveints and sendints. You need this parameter when
+ | calling these routines. Note that for many calls I can use
+ | the variable 'smallidx' which is exactly the number of bits, and
+ | So I don't need to call 'sizeofints for those calls.
+*/
+
+static int sizeofints( const int num_of_ints, unsigned int sizes[]) {
+ int i, num;
+ unsigned int num_of_bytes, num_of_bits, bytes[32], bytecnt, tmp;
+ num_of_bytes = 1;
+ bytes[0] = 1;
+ num_of_bits = 0;
+ for (i=0; i < num_of_ints; i++) {
+ tmp = 0;
+ for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++) {
+ tmp = bytes[bytecnt] * sizes[i] + tmp;
+ bytes[bytecnt] = tmp & 0xff;
+ tmp >>= 8;
+ }
+ while (tmp != 0) {
+ bytes[bytecnt++] = tmp & 0xff;
+ tmp >>= 8;
+ }
+ num_of_bytes = bytecnt;
+ }
+ num = 1;
+ num_of_bytes--;
+ while (bytes[num_of_bytes] >= num) {
+ num_of_bits++;
+ num *= 2;
+ }
+ return num_of_bits + num_of_bytes * 8;
+
+}
+
+/*____________________________________________________________________________
+ |
+ | sendints - send a small set of small integers in compressed format
+ |
+ | this routine is used internally by xdr3dfcoord, to send a set of
+ | small integers to the buffer.
+ | Multiplication with fixed (specified maximum ) sizes is used to get
+ | to one big, multibyte integer. Allthough the routine could be
+ | modified to handle sizes bigger than 16777216, or more than just
+ | a few integers, this is not done, because the gain in compression
+ | isn't worth the effort. Note that overflowing the multiplication
+ | or the byte buffer (32 bytes) is unchecked and causes bad results.
+ |
+ */
+
+static void sendints(int buf[], const int num_of_ints, const int num_of_bits,
+ unsigned int sizes[], unsigned int nums[]) {
+
+ int i;
+ unsigned int bytes[32], num_of_bytes, bytecnt, tmp;
+
+ tmp = nums[0];
+ num_of_bytes = 0;
+ do {
+ bytes[num_of_bytes++] = tmp & 0xff;
+ tmp >>= 8;
+ } while (tmp != 0);
+
+ for (i = 1; i < num_of_ints; i++) {
+ if (nums[i] >= sizes[i]) {
+ fprintf(stderr,"major breakdown in sendints num %d doesn't "
+ "match size %d\n", nums[i], sizes[i]);
+ exit(1);
+ }
+ /* use one step multiply */
+ tmp = nums[i];
+ for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++) {
+ tmp = bytes[bytecnt] * sizes[i] + tmp;
+ bytes[bytecnt] = tmp & 0xff;
+ tmp >>= 8;
+ }
+ while (tmp != 0) {
+ bytes[bytecnt++] = tmp & 0xff;
+ tmp >>= 8;
+ }
+ num_of_bytes = bytecnt;
+ }
+ if (num_of_bits >= num_of_bytes * 8) {
+ for (i = 0; i < num_of_bytes; i++) {
+ sendbits(buf, 8, bytes[i]);
+ }
+ sendbits(buf, num_of_bits - num_of_bytes * 8, 0);
+ } else {
+ for (i = 0; i < num_of_bytes-1; i++) {
+ sendbits(buf, 8, bytes[i]);
+ }
+ sendbits(buf, num_of_bits- (num_of_bytes -1) * 8, bytes[i]);
+ }
+}
+
+
+/*___________________________________________________________________________
+ |
+ | receivebits - decode number from buf using specified number of bits
+ |
+ | extract the number of bits from the array buf and construct an integer
+ | from it. Return that value.
+ |
+*/
+
+static int receivebits(int buf[], int num_of_bits) {
+
+ int cnt, num;
+ unsigned int lastbits, lastbyte;
+ unsigned char * cbuf;
+ int mask = (1 << num_of_bits) -1;
+
+ cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
+ cnt = buf[0];
+ lastbits = (unsigned int) buf[1];
+ lastbyte = (unsigned int) buf[2];
+
+ num = 0;
+ while (num_of_bits >= 8) {
+ lastbyte = ( lastbyte << 8 ) | cbuf[cnt++];
+ num |= (lastbyte >> lastbits) << (num_of_bits - 8);
+ num_of_bits -=8;
+ }
+ if (num_of_bits > 0) {
+ if (lastbits < num_of_bits) {
+ lastbits += 8;
+ lastbyte = (lastbyte << 8) | cbuf[cnt++];
+ }
+ lastbits -= num_of_bits;
+ num |= (lastbyte >> lastbits) & ((1 << num_of_bits) -1);
+ }
+ num &= mask;
+ buf[0] = cnt;
+ buf[1] = lastbits;
+ buf[2] = lastbyte;
+ return num;
+}
+
+/*____________________________________________________________________________
+ |
+ | receiveints - decode 'small' integers from the buf array
+ |
+ | this routine is the inverse from sendints() and decodes the small integers
+ | written to buf by calculating the remainder and doing divisions with
+ | the given sizes[]. You need to specify the total number of bits to be
+ | used from buf in num_of_bits.
+ |
+*/
+
+static void receiveints(int buf[], const int num_of_ints, int num_of_bits,
+ unsigned int sizes[], int nums[]) {
+ int bytes[32];
+ int i, j, num_of_bytes, p, num;
+
+ bytes[1] = bytes[2] = bytes[3] = 0;
+ num_of_bytes = 0;
+ while (num_of_bits > 8) {
+ bytes[num_of_bytes++] = receivebits(buf, 8);
+ num_of_bits -= 8;
+ }
+ if (num_of_bits > 0) {
+ bytes[num_of_bytes++] = receivebits(buf, num_of_bits);
+ }
+ for (i = num_of_ints-1; i > 0; i--) {
+ num = 0;
+ for (j = num_of_bytes-1; j >=0; j--) {
+ num = (num << 8) | bytes[j];
+ p = num / sizes[i];
+ bytes[j] = p;
+ num = num - p * sizes[i];
+ }
+ nums[i] = num;
+ }
+ nums[0] = bytes[0] | (bytes[1] << 8) | (bytes[2] << 16) | (bytes[3] << 24);
+}
+
+/*____________________________________________________________________________
+ |
+ | xdr3dfcoord - read or write compressed 3d coordinates to xdr file.
+ |
+ | this routine reads or writes (depending on how you opened the file with
+ | xdropen() ) a large number of 3d coordinates (stored in *fp).
+ | The number of coordinates triplets to write is given by *size. On
+ | read this number may be zero, in which case it reads as many as were written
+ | or it may specify the number if triplets to read (which should match the
+ | number written).
+ | Compression is achieved by first converting all floating numbers to integer
+ | using multiplication by *precision and rounding to the nearest integer.
+ | Then the minimum and maximum value are calculated to determine the range.
+ | The limited range of integers so found, is used to compress the coordinates.
+ | In addition the differences between succesive coordinates is calculated.
+ | If the difference happens to be 'small' then only the difference is saved,
+ | compressing the data even more. The notion of 'small' is changed dynamically
+ | and is enlarged or reduced whenever needed or possible.
+ | Extra compression is achieved in the case of GROMOS and coordinates of
+ | water molecules. GROMOS first writes out the Oxygen position, followed by
+ | the two hydrogens. In order to make the differences smaller (and thereby
+ | compression the data better) the order is changed into first one hydrogen
+ | then the oxygen, followed by the other hydrogen. This is rather special, but
+ | it shouldn't harm in the general case.
+ |
+ */
+
+int xdr3dfcoord(XDR *xdrs, float *fp, int *size, float *precision) {
+
+
+ static int *ip = NULL;
+ static int oldsize;
+ static int *buf;
+
+ int minint[3], maxint[3], mindiff, *lip, diff;
+ int lint1, lint2, lint3, oldlint1, oldlint2, oldlint3, smallidx;
+ int minidx, maxidx;
+ unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3, *luip;
+ int flag, k;
+ int small, smaller, larger, i, is_small, is_smaller, run, prevrun;
+ float *lfp, lf;
+ int tmp, *thiscoord, prevcoord[3];
+ unsigned int tmpcoord[30];
+
+ int bufsize, xdrid, lsize;
+ unsigned int bitsize;
+ float inv_precision;
+ int errval = 1;
+
+ /* find out if xdrs is opened for reading or for writing */
+ xdrid = 0;
+ while (xdridptr[xdrid] != xdrs) {
+ xdrid++;
+ if (xdrid >= MAXID) {
+ fprintf(stderr, "xdr error. no open xdr stream\n");
+ exit (1);
+ }
+ }
+ if (xdrmodes[xdrid] == 'w') {
+
+ /* xdrs is open for writing */
+
+ if (xdr_int(xdrs, size) == 0)
+ return 0;
+ size3 = *size * 3;
+ /* when the number of coordinates is small, don't try to compress; just
+ * write them as floats using xdr_vector
+ */
+ if (*size <= 9 ) {
+ return (xdr_vector(xdrs, (char *) fp, size3, sizeof(*fp),
+ (xdrproc_t)xdr_float));
+ }
+
+ xdr_float(xdrs, precision);
+ if (ip == NULL) {
+ ip = (int *)malloc(size3 * sizeof(*ip));
+ if (ip == NULL) {
+ fprintf(stderr,"malloc failed\n");
+ exit(1);
+ }
+ bufsize = size3 * 1.2;
+ buf = (int *)malloc(bufsize * sizeof(*buf));
+ if (buf == NULL) {
+ fprintf(stderr,"malloc failed\n");
+ exit(1);
+ }
+ oldsize = *size;
+ } else if (*size > oldsize) {
+ ip = (int *)realloc(ip, size3 * sizeof(*ip));
+ if (ip == NULL) {
+ fprintf(stderr,"malloc failed\n");
+ exit(1);
+ }
+ bufsize = size3 * 1.2;
+ buf = (int *)realloc(buf, bufsize * sizeof(*buf));
+ if (buf == NULL) {
+ fprintf(stderr,"malloc failed\n");
+ exit(1);
+ }
+ oldsize = *size;
+ }
+ /* buf[0-2] are special and do not contain actual data */
+ buf[0] = buf[1] = buf[2] = 0;
+ minint[0] = minint[1] = minint[2] = INT_MAX;
+ maxint[0] = maxint[1] = maxint[2] = INT_MIN;
+ prevrun = -1;
+ lfp = fp;
+ lip = ip;
+ mindiff = INT_MAX;
+ oldlint1 = oldlint2 = oldlint3 = 0;
+ while(lfp < fp + size3 ) {
+ /* find nearest integer */
+ if (*lfp >= 0.0)
+ lf = *lfp * *precision + 0.5;
+ else
+ lf = *lfp * *precision - 0.5;
+ if (fabs(lf) > MAXABS) {
+ /* scaling would cause overflow */
+ errval = 0;
+ }
+ lint1 = lf;
+ if (lint1 < minint[0]) minint[0] = lint1;
+ if (lint1 > maxint[0]) maxint[0] = lint1;
+ *lip++ = lint1;
+ lfp++;
+ if (*lfp >= 0.0)
+ lf = *lfp * *precision + 0.5;
+ else
+ lf = *lfp * *precision - 0.5;
+ if (fabs(lf) > MAXABS) {
+ /* scaling would cause overflow */
+ errval = 0;
+ }
+ lint2 = lf;
+ if (lint2 < minint[1]) minint[1] = lint2;
+ if (lint2 > maxint[1]) maxint[1] = lint2;
+ *lip++ = lint2;
+ lfp++;
+ if (*lfp >= 0.0)
+ lf = *lfp * *precision + 0.5;
+ else
+ lf = *lfp * *precision - 0.5;
+ if (fabs(lf) > MAXABS) {
+ /* scaling would cause overflow */
+ errval = 0;
+ }
+ lint3 = lf;
+ if (lint3 < minint[2]) minint[2] = lint3;
+ if (lint3 > maxint[2]) maxint[2] = lint3;
+ *lip++ = lint3;
+ lfp++;
+ diff = abs(oldlint1-lint1)+abs(oldlint2-lint2)+abs(oldlint3-lint3);
+ if (diff < mindiff && lfp > fp + 3)
+ mindiff = diff;
+ oldlint1 = lint1;
+ oldlint2 = lint2;
+ oldlint3 = lint3;
+ }
+ xdr_int(xdrs, &(minint[0]));
+ xdr_int(xdrs, &(minint[1]));
+ xdr_int(xdrs, &(minint[2]));
+
+ xdr_int(xdrs, &(maxint[0]));
+ xdr_int(xdrs, &(maxint[1]));
+ xdr_int(xdrs, &(maxint[2]));
+
+ if ((float)maxint[0] - (float)minint[0] >= MAXABS ||
+ (float)maxint[1] - (float)minint[1] >= MAXABS ||
+ (float)maxint[2] - (float)minint[2] >= MAXABS) {
+ /* turning value in unsigned by subtracting minint
+ * would cause overflow
+ */
+ errval = 0;
+ }
+ sizeint[0] = maxint[0] - minint[0]+1;
+ sizeint[1] = maxint[1] - minint[1]+1;
+ sizeint[2] = maxint[2] - minint[2]+1;
+
+ /* check if one of the sizes is to big to be multiplied */
+ if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) {
+ bitsizeint[0] = sizeofint(sizeint[0]);
+ bitsizeint[1] = sizeofint(sizeint[1]);
+ bitsizeint[2] = sizeofint(sizeint[2]);
+ bitsize = 0; /* flag the use of large sizes */
+ } else {
+ bitsize = sizeofints(3, sizeint);
+ }
+ lip = ip;
+ luip = (unsigned int *) ip;
+ smallidx = FIRSTIDX;
+ while (smallidx < LASTIDX && magicints[smallidx] < mindiff) {
+ smallidx++;
+ }
+ xdr_int(xdrs, &smallidx);
+ maxidx = MIN(LASTIDX, smallidx + 8) ;
+ minidx = maxidx - 8; /* often this equal smallidx */
+ smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 2;
+ small = magicints[smallidx] / 2;
+ sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
+ larger = magicints[maxidx] / 2;
+ i = 0;
+ while (i < *size) {
+ is_small = 0;
+ thiscoord = (int *)(luip) + i * 3;
+ if (smallidx < maxidx && i >= 1 &&
+ abs(thiscoord[0] - prevcoord[0]) < larger &&
+ abs(thiscoord[1] - prevcoord[1]) < larger &&
+ abs(thiscoord[2] - prevcoord[2]) < larger) {
+ is_smaller = 1;
+ } else if (smallidx > minidx) {
+ is_smaller = -1;
+ } else {
+ is_smaller = 0;
+ }
+ if (i + 1 < *size) {
+ if (abs(thiscoord[0] - thiscoord[3]) < small &&
+ abs(thiscoord[1] - thiscoord[4]) < small &&
+ abs(thiscoord[2] - thiscoord[5]) < small) {
+ /* interchange first with second atom for better
+ * compression of water molecules
+ */
+ tmp = thiscoord[0]; thiscoord[0] = thiscoord[3];
+ thiscoord[3] = tmp;
+ tmp = thiscoord[1]; thiscoord[1] = thiscoord[4];
+ thiscoord[4] = tmp;
+ tmp = thiscoord[2]; thiscoord[2] = thiscoord[5];
+ thiscoord[5] = tmp;
+ is_small = 1;
+ }
+
+ }
+ tmpcoord[0] = thiscoord[0] - minint[0];
+ tmpcoord[1] = thiscoord[1] - minint[1];
+ tmpcoord[2] = thiscoord[2] - minint[2];
+ if (bitsize == 0) {
+ sendbits(buf, bitsizeint[0], tmpcoord[0]);
+ sendbits(buf, bitsizeint[1], tmpcoord[1]);
+ sendbits(buf, bitsizeint[2], tmpcoord[2]);
+ } else {
+ sendints(buf, 3, bitsize, sizeint, tmpcoord);
+ }
+ prevcoord[0] = thiscoord[0];
+ prevcoord[1] = thiscoord[1];
+ prevcoord[2] = thiscoord[2];
+ thiscoord = thiscoord + 3;
+ i++;
+
+ run = 0;
+ if (is_small == 0 && is_smaller == -1)
+ is_smaller = 0;
+ while (is_small && run < 8*3) {
+ if (is_smaller == -1 && (
+ SQR(thiscoord[0] - prevcoord[0]) +
+ SQR(thiscoord[1] - prevcoord[1]) +
+ SQR(thiscoord[2] - prevcoord[2]) >= smaller * smaller)) {
+ is_smaller = 0;
+ }
+
+ tmpcoord[run++] = thiscoord[0] - prevcoord[0] + small;
+ tmpcoord[run++] = thiscoord[1] - prevcoord[1] + small;
+ tmpcoord[run++] = thiscoord[2] - prevcoord[2] + small;
+
+ prevcoord[0] = thiscoord[0];
+ prevcoord[1] = thiscoord[1];
+ prevcoord[2] = thiscoord[2];
+
+ i++;
+ thiscoord = thiscoord + 3;
+ is_small = 0;
+ if (i < *size &&
+ abs(thiscoord[0] - prevcoord[0]) < small &&
+ abs(thiscoord[1] - prevcoord[1]) < small &&
+ abs(thiscoord[2] - prevcoord[2]) < small) {
+ is_small = 1;
+ }
+ }
+ if (run != prevrun || is_smaller != 0) {
+ prevrun = run;
+ sendbits(buf, 1, 1); /* flag the change in run-length */
+ sendbits(buf, 5, run+is_smaller+1);
+ } else {
+ sendbits(buf, 1, 0); /* flag the fact that runlength did not change */
+ }
+ for (k=0; k < run; k+=3) {
+ sendints(buf, 3, smallidx, sizesmall, &tmpcoord[k]);
+ }
+ if (is_smaller != 0) {
+ smallidx += is_smaller;
+ if (is_smaller < 0) {
+ small = smaller;
+ smaller = magicints[smallidx-1] / 2;
+ } else {
+ smaller = small;
+ small = magicints[smallidx] / 2;
+ }
+ sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
+ }
+ }
+ if (buf[1] != 0) buf[0]++;;
+ xdr_int(xdrs, &(buf[0])); /* buf[0] holds the length in bytes */
+ return errval * (xdr_opaque(xdrs, (caddr_t)&(buf[3]), (u_int)buf[0]));
+ } else {
+
+ /* xdrs is open for reading */
+
+ if (xdr_int(xdrs, &lsize) == 0)
+ return 0;
+ if (*size != 0 && lsize != *size) {
+ fprintf(stderr, "wrong number of coordinates in xdr3dfcoor; "
+ "%d arg vs %d in file", *size, lsize);
+ }
+ *size = lsize;
+ size3 = *size * 3;
+ if (*size <= 9) {
+ return (xdr_vector(xdrs, (char *) fp, size3, sizeof(*fp),
+ (xdrproc_t)xdr_float));
+ }
+ xdr_float(xdrs, precision);
+ if (ip == NULL) {
+ ip = (int *)malloc(size3 * sizeof(*ip));
+ if (ip == NULL) {
+ fprintf(stderr,"malloc failed\n");
+ exit(1);
+ }
+ bufsize = size3 * 1.2;
+ buf = (int *)malloc(bufsize * sizeof(*buf));
+ if (buf == NULL) {
+ fprintf(stderr,"malloc failed\n");
+ exit(1);
+ }
+ oldsize = *size;
+ } else if (*size > oldsize) {
+ ip = (int *)realloc(ip, size3 * sizeof(*ip));
+ if (ip == NULL) {
+ fprintf(stderr,"malloc failed\n");
+ exit(1);
+ }
+ bufsize = size3 * 1.2;
+ buf = (int *)realloc(buf, bufsize * sizeof(*buf));
+ if (buf == NULL) {
+ fprintf(stderr,"malloc failed\n");
+ exit(1);
+ }
+ oldsize = *size;
+ }
+ buf[0] = buf[1] = buf[2] = 0;
+
+ xdr_int(xdrs, &(minint[0]));
+ xdr_int(xdrs, &(minint[1]));
+ xdr_int(xdrs, &(minint[2]));
+
+ xdr_int(xdrs, &(maxint[0]));
+ xdr_int(xdrs, &(maxint[1]));
+ xdr_int(xdrs, &(maxint[2]));
+
+ sizeint[0] = maxint[0] - minint[0]+1;
+ sizeint[1] = maxint[1] - minint[1]+1;
+ sizeint[2] = maxint[2] - minint[2]+1;
+
+ /* check if one of the sizes is to big to be multiplied */
+ if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) {
+ bitsizeint[0] = sizeofint(sizeint[0]);
+ bitsizeint[1] = sizeofint(sizeint[1]);
+ bitsizeint[2] = sizeofint(sizeint[2]);
+ bitsize = 0; /* flag the use of large sizes */
+ } else {
+ bitsize = sizeofints(3, sizeint);
+ }
+
+ xdr_int(xdrs, &smallidx);
+ maxidx = MIN(LASTIDX, smallidx + 8) ;
+ minidx = maxidx - 8; /* often this equal smallidx */
+ smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 2;
+ small = magicints[smallidx] / 2;
+ sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ;
+ larger = magicints[maxidx];
+
+ /* buf[0] holds the length in bytes */
+
+ if (xdr_int(xdrs, &(buf[0])) == 0)
+ return 0;
+ if (xdr_opaque(xdrs, (caddr_t)&(buf[3]), (u_int)buf[0]) == 0)
+ return 0;
+ buf[0] = buf[1] = buf[2] = 0;
+
+ lfp = fp;
+ inv_precision = 1.0 / * precision;
+ run = 0;
+ i = 0;
+ lip = ip;
+ while ( i < lsize ) {
+ thiscoord = (int *)(lip) + i * 3;
+
+ if (bitsize == 0) {
+ thiscoord[0] = receivebits(buf, bitsizeint[0]);
+ thiscoord[1] = receivebits(buf, bitsizeint[1]);
+ thiscoord[2] = receivebits(buf, bitsizeint[2]);
+ } else {
+ receiveints(buf, 3, bitsize, sizeint, thiscoord);
+ }
+
+ i++;
+ thiscoord[0] += minint[0];
+ thiscoord[1] += minint[1];
+ thiscoord[2] += minint[2];
+
+ prevcoord[0] = thiscoord[0];
+ prevcoord[1] = thiscoord[1];
+ prevcoord[2] = thiscoord[2];
+
+
+ flag = receivebits(buf, 1);
+ is_smaller = 0;
+ if (flag == 1) {
+ run = receivebits(buf, 5);
+ is_smaller = run % 3;
+ run -= is_smaller;
+ is_smaller--;
+ }
+ if (run > 0) {
+ thiscoord += 3;
+ for (k = 0; k < run; k+=3) {
+ receiveints(buf, 3, smallidx, sizesmall, thiscoord);
+ i++;
+ thiscoord[0] += prevcoord[0] - small;
+ thiscoord[1] += prevcoord[1] - small;
+ thiscoord[2] += prevcoord[2] - small;
+ if (k == 0) {
+ /* interchange first with second atom for better
+ * compression of water molecules
+ */
+ tmp = thiscoord[0]; thiscoord[0] = prevcoord[0];
+ prevcoord[0] = tmp;
+ tmp = thiscoord[1]; thiscoord[1] = prevcoord[1];
+ prevcoord[1] = tmp;
+ tmp = thiscoord[2]; thiscoord[2] = prevcoord[2];
+ prevcoord[2] = tmp;
+ *lfp++ = prevcoord[0] * inv_precision;
+ *lfp++ = prevcoord[1] * inv_precision;
+ *lfp++ = prevcoord[2] * inv_precision;
+ } else {
+ prevcoord[0] = thiscoord[0];
+ prevcoord[1] = thiscoord[1];
+ prevcoord[2] = thiscoord[2];
+ }
+ *lfp++ = thiscoord[0] * inv_precision;
+ *lfp++ = thiscoord[1] * inv_precision;
+ *lfp++ = thiscoord[2] * inv_precision;
+ }
+ } else {
+ *lfp++ = thiscoord[0] * inv_precision;
+ *lfp++ = thiscoord[1] * inv_precision;
+ *lfp++ = thiscoord[2] * inv_precision;
+ }
+ smallidx += is_smaller;
+ if (is_smaller < 0) {
+ small = smaller;
+ if (smallidx > FIRSTIDX) {
+ smaller = magicints[smallidx - 1] /2;
+ } else {
+ smaller = 0;
+ }
+ } else if (is_smaller > 0) {
+ smaller = small;
+ small = magicints[smallidx] / 2;
+ }
+ sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ;
+ }
+ }
+ return 1;
+}
+
+
+
--- /dev/null
+/*
+ * Sun RPC is a product of Sun Microsystems, Inc. and is provided for
+ * unrestricted use provided that this legend is included on all tape
+ * media and as a part of the software program in whole or part. Users
+ * may copy or modify Sun RPC without charge, but are not authorized
+ * to license or distribute it to anyone else except as part of a product or
+ * program developed by the user.
+ *
+ * SUN RPC IS PROVIDED AS IS WITH NO WARRANTIES OF ANY KIND INCLUDING THE
+ * WARRANTIES OF DESIGN, MERCHANTIBILITY AND FITNESS FOR A PARTICULAR
+ * PURPOSE, OR ARISING FROM A COURSE OF DEALING, USAGE OR TRADE PRACTICE.
+ *
+ * Sun RPC is provided with no support and without any obligation on the
+ * part of Sun Microsystems, Inc. to assist in its use, correction,
+ * modification or enhancement.
+ *
+ * SUN MICROSYSTEMS, INC. SHALL HAVE NO LIABILITY WITH RESPECT TO THE
+ * INFRINGEMENT OF COPYRIGHTS, TRADE SECRETS OR ANY PATENTS BY SUN RPC
+ * OR ANY PART THEREOF.
+ *
+ * In no event will Sun Microsystems, Inc. be liable for any lost revenue
+ * or profits or other special, indirect and consequential damages, even if
+ * Sun has been advised of the possibility of such damages.
+ *
+ * Sun Microsystems, Inc.
+ * 2550 Garcia Avenue
+ * Mountain View, California 94043
+ */
+/* fixincludes should not add extern "C" to this file */
+/*
+ * Rpc additions to <sys/types.h>
+ */
+#ifndef _RPC_TYPES_H
+#define _RPC_TYPES_H 1
+
+typedef int bool_t;
+typedef int enum_t;
+/* This needs to be changed to uint32_t in the future */
+typedef unsigned long rpcprog_t;
+typedef unsigned long rpcvers_t;
+typedef unsigned long rpcproc_t;
+typedef unsigned long rpcprot_t;
+typedef unsigned long rpcport_t;
+
+#define __dontcare__ -1
+
+#ifndef FALSE
+# define FALSE (0)
+#endif
+
+#ifndef TRUE
+# define TRUE (1)
+#endif
+
+#ifndef NULL
+# define NULL 0
+#endif
+
+#include <stdlib.h> /* For malloc decl. */
+#define mem_alloc(bsize) malloc(bsize)
+/*
+ * XXX: This must not use the second argument, or code in xdr_array.c needs
+ * to be modified.
+ */
+#define mem_free(ptr, bsize) free(ptr)
+
+#ifndef makedev /* ie, we haven't already included it */
+#include <sys/types.h>
+#endif
+
+#ifndef __u_char_defined
+typedef __u_char u_char;
+typedef __u_short u_short;
+typedef __u_int u_int;
+typedef __u_long u_long;
+typedef __quad_t quad_t;
+typedef __u_quad_t u_quad_t;
+typedef __fsid_t fsid_t;
+# define __u_char_defined
+#endif
+#ifndef __daddr_t_defined
+typedef __daddr_t daddr_t;
+typedef __caddr_t caddr_t;
+# define __daddr_t_defined
+#endif
+
+#include <sys/time.h>
+#include <sys/param.h>
+
+#include <netinet/in.h>
+
+#ifndef INADDR_LOOPBACK
+#define INADDR_LOOPBACK (u_long)0x7F000001
+#endif
+#ifndef MAXHOSTNAMELEN
+#define MAXHOSTNAMELEN 64
+#endif
+
+#endif /* rpc/types.h */
--- /dev/null
+divert(-1)
+undefine(`len')
+#
+# append an underscore to FORTRAN function names
+#
+define(`FUNCTION',`$1_')
+#
+# FORTRAN character strings are passed as follows:
+# a pointer to the base of the string is passed in the normal
+# argument list, and the length is passed by value as an extra
+# argument, after all of the other arguments.
+#
+define(`ARGS',`($1`'undivert(1))')
+define(`SAVE',`divert(1)$1`'divert(0)')
+define(`STRING_ARG',`$1_ptr`'SAVE(`, $1_len')')
+define(`STRING_ARG_DECL',`char * $1_ptr; int $1_len')
+define(`STRING_LEN',`$1_len')
+define(`STRING_PTR',`$1_ptr')
+divert(0)
--- /dev/null
+# define INTUSE(name) name
+# define INTDEF(name)
+/* @(#)xdr.c 2.1 88/07/29 4.0 RPCSRC */
+/*
+ * Sun RPC is a product of Sun Microsystems, Inc. and is provided for
+ * unrestricted use provided that this legend is included on all tape
+ * media and as a part of the software program in whole or part. Users
+ * may copy or modify Sun RPC without charge, but are not authorized
+ * to license or distribute it to anyone else except as part of a product or
+ * program developed by the user.
+ *
+ * SUN RPC IS PROVIDED AS IS WITH NO WARRANTIES OF ANY KIND INCLUDING THE
+ * WARRANTIES OF DESIGN, MERCHANTIBILITY AND FITNESS FOR A PARTICULAR
+ * PURPOSE, OR ARISING FROM A COURSE OF DEALING, USAGE OR TRADE PRACTICE.
+ *
+ * Sun RPC is provided with no support and without any obligation on the
+ * part of Sun Microsystems, Inc. to assist in its use, correction,
+ * modification or enhancement.
+ *
+ * SUN MICROSYSTEMS, INC. SHALL HAVE NO LIABILITY WITH RESPECT TO THE
+ * INFRINGEMENT OF COPYRIGHTS, TRADE SECRETS OR ANY PATENTS BY SUN RPC
+ * OR ANY PART THEREOF.
+ *
+ * In no event will Sun Microsystems, Inc. be liable for any lost revenue
+ * or profits or other special, indirect and consequential damages, even if
+ * Sun has been advised of the possibility of such damages.
+ *
+ * Sun Microsystems, Inc.
+ * 2550 Garcia Avenue
+ * Mountain View, California 94043
+ */
+#if !defined(lint) && defined(SCCSIDS)
+static char sccsid[] = "@(#)xdr.c 1.35 87/08/12";
+#endif
+
+/*
+ * xdr.c, Generic XDR routines implementation.
+ *
+ * Copyright (C) 1986, Sun Microsystems, Inc.
+ *
+ * These are the "generic" xdr routines used to serialize and de-serialize
+ * most common data items. See xdr.h for more info on the interface to
+ * xdr.
+ */
+
+#include <stdio.h>
+#include <limits.h>
+#include <string.h>
+#include <libintl.h>
+
+#include "types.h"
+#include "xdr.h"
+
+#ifdef USE_IN_LIBIO
+# include <wchar.h>
+#endif
+
+/*
+ * constants specific to the xdr "protocol"
+ */
+#define XDR_FALSE ((long) 0)
+#define XDR_TRUE ((long) 1)
+#define LASTUNSIGNED ((u_int) 0-1)
+
+/*
+ * for unit alignment
+ */
+static const char xdr_zero[BYTES_PER_XDR_UNIT] = {0, 0, 0, 0};
+
+/*
+ * Free a data structure using XDR
+ * Not a filter, but a convenient utility nonetheless
+ */
+void
+xdr_free (xdrproc_t proc, char *objp)
+{
+ XDR x;
+
+ x.x_op = XDR_FREE;
+ (*proc) (&x, objp);
+}
+
+/*
+ * XDR nothing
+ */
+bool_t
+xdr_void (void)
+{
+ return TRUE;
+}
+INTDEF(xdr_void)
+
+/*
+ * XDR integers
+ */
+bool_t
+xdr_int (XDR *xdrs, int *ip)
+{
+
+#if INT_MAX < LONG_MAX
+ long l;
+
+ switch (xdrs->x_op)
+ {
+ case XDR_ENCODE:
+ l = (long) *ip;
+ return XDR_PUTLONG (xdrs, &l);
+
+ case XDR_DECODE:
+ if (!XDR_GETLONG (xdrs, &l))
+ {
+ return FALSE;
+ }
+ *ip = (int) l;
+ case XDR_FREE:
+ return TRUE;
+ }
+ return FALSE;
+#elif INT_MAX == LONG_MAX
+ return INTUSE(xdr_long) (xdrs, (long *) ip);
+#elif INT_MAX == SHRT_MAX
+ return INTUSE(xdr_short) (xdrs, (short *) ip);
+#else
+#error unexpected integer sizes in_xdr_int()
+#endif
+}
+INTDEF(xdr_int)
+
+/*
+ * XDR unsigned integers
+ */
+bool_t
+xdr_u_int (XDR *xdrs, u_int *up)
+{
+#if UINT_MAX < ULONG_MAX
+ long l;
+
+ switch (xdrs->x_op)
+ {
+ case XDR_ENCODE:
+ l = (u_long) * up;
+ return XDR_PUTLONG (xdrs, &l);
+
+ case XDR_DECODE:
+ if (!XDR_GETLONG (xdrs, &l))
+ {
+ return FALSE;
+ }
+ *up = (u_int) (u_long) l;
+ case XDR_FREE:
+ return TRUE;
+ }
+ return FALSE;
+#elif UINT_MAX == ULONG_MAX
+ return INTUSE(xdr_u_long) (xdrs, (u_long *) up);
+#elif UINT_MAX == USHRT_MAX
+ return INTUSE(xdr_short) (xdrs, (short *) up);
+#else
+#error unexpected integer sizes in_xdr_u_int()
+#endif
+}
+INTDEF(xdr_u_int)
+
+/*
+ * XDR long integers
+ * The definition of xdr_long() is kept for backward
+ * compatibility. Instead xdr_int() should be used.
+ */
+bool_t
+xdr_long (XDR *xdrs, long *lp)
+{
+
+ if (xdrs->x_op == XDR_ENCODE
+ && (sizeof (int32_t) == sizeof (long)
+ || (int32_t) *lp == *lp))
+ return XDR_PUTLONG (xdrs, lp);
+
+ if (xdrs->x_op == XDR_DECODE)
+ return XDR_GETLONG (xdrs, lp);
+
+ if (xdrs->x_op == XDR_FREE)
+ return TRUE;
+
+ return FALSE;
+}
+INTDEF(xdr_long)
+
+/*
+ * XDR unsigned long integers
+ * The definition of xdr_u_long() is kept for backward
+ * compatibility. Instead xdr_u_int() should be used.
+ */
+bool_t
+xdr_u_long (XDR *xdrs, u_long *ulp)
+{
+ switch (xdrs->x_op)
+ {
+ case XDR_DECODE:
+ {
+ long int tmp;
+
+ if (XDR_GETLONG (xdrs, &tmp) == FALSE)
+ return FALSE;
+
+ *ulp = (uint32_t) tmp;
+ return TRUE;
+ }
+
+ case XDR_ENCODE:
+ if (sizeof (uint32_t) != sizeof (u_long)
+ && (uint32_t) *ulp != *ulp)
+ return FALSE;
+
+ return XDR_PUTLONG (xdrs, (long *) ulp);
+
+ case XDR_FREE:
+ return TRUE;
+ }
+ return FALSE;
+}
+INTDEF(xdr_u_long)
+
+/*
+ * XDR hyper integers
+ * same as xdr_u_hyper - open coded to save a proc call!
+ */
+bool_t
+xdr_hyper (XDR *xdrs, quad_t *llp)
+{
+ long int t1, t2;
+
+ if (xdrs->x_op == XDR_ENCODE)
+ {
+ t1 = (long) ((*llp) >> 32);
+ t2 = (long) (*llp);
+ return (XDR_PUTLONG(xdrs, &t1) && XDR_PUTLONG(xdrs, &t2));
+ }
+
+ if (xdrs->x_op == XDR_DECODE)
+ {
+ if (!XDR_GETLONG(xdrs, &t1) || !XDR_GETLONG(xdrs, &t2))
+ return FALSE;
+ *llp = ((quad_t) t1) << 32;
+ *llp |= (uint32_t) t2;
+ return TRUE;
+ }
+
+ if (xdrs->x_op == XDR_FREE)
+ return TRUE;
+
+ return FALSE;
+}
+INTDEF(xdr_hyper)
+
+
+/*
+ * XDR hyper integers
+ * same as xdr_hyper - open coded to save a proc call!
+ */
+bool_t
+xdr_u_hyper (XDR *xdrs, u_quad_t *ullp)
+{
+ long int t1, t2;
+
+ if (xdrs->x_op == XDR_ENCODE)
+ {
+ t1 = (unsigned long) ((*ullp) >> 32);
+ t2 = (unsigned long) (*ullp);
+ return (XDR_PUTLONG(xdrs, &t1) && XDR_PUTLONG(xdrs, &t2));
+ }
+
+ if (xdrs->x_op == XDR_DECODE)
+ {
+ if (!XDR_GETLONG(xdrs, &t1) || !XDR_GETLONG(xdrs, &t2))
+ return FALSE;
+ *ullp = ((u_quad_t) t1) << 32;
+ *ullp |= (uint32_t) t2;
+ return TRUE;
+ }
+
+ if (xdrs->x_op == XDR_FREE)
+ return TRUE;
+
+ return FALSE;
+}
+INTDEF(xdr_u_hyper)
+
+bool_t
+xdr_longlong_t (XDR *xdrs, quad_t *llp)
+{
+ return INTUSE(xdr_hyper) (xdrs, llp);
+}
+
+bool_t
+xdr_u_longlong_t (XDR *xdrs, u_quad_t *ullp)
+{
+ return INTUSE(xdr_u_hyper) (xdrs, ullp);
+}
+
+/*
+ * XDR short integers
+ */
+bool_t
+xdr_short (XDR *xdrs, short *sp)
+{
+ long l;
+
+ switch (xdrs->x_op)
+ {
+ case XDR_ENCODE:
+ l = (long) *sp;
+ return XDR_PUTLONG (xdrs, &l);
+
+ case XDR_DECODE:
+ if (!XDR_GETLONG (xdrs, &l))
+ {
+ return FALSE;
+ }
+ *sp = (short) l;
+ return TRUE;
+
+ case XDR_FREE:
+ return TRUE;
+ }
+ return FALSE;
+}
+INTDEF(xdr_short)
+
+/*
+ * XDR unsigned short integers
+ */
+bool_t
+xdr_u_short (XDR *xdrs, u_short *usp)
+{
+ long l;
+
+ switch (xdrs->x_op)
+ {
+ case XDR_ENCODE:
+ l = (u_long) * usp;
+ return XDR_PUTLONG (xdrs, &l);
+
+ case XDR_DECODE:
+ if (!XDR_GETLONG (xdrs, &l))
+ {
+ return FALSE;
+ }
+ *usp = (u_short) (u_long) l;
+ return TRUE;
+
+ case XDR_FREE:
+ return TRUE;
+ }
+ return FALSE;
+}
+INTDEF(xdr_u_short)
+
+
+/*
+ * XDR a char
+ */
+bool_t
+xdr_char (XDR *xdrs, char *cp)
+{
+ int i;
+
+ i = (*cp);
+ if (!INTUSE(xdr_int) (xdrs, &i))
+ {
+ return FALSE;
+ }
+ *cp = i;
+ return TRUE;
+}
+
+/*
+ * XDR an unsigned char
+ */
+bool_t
+xdr_u_char (XDR *xdrs, u_char *cp)
+{
+ u_int u;
+
+ u = (*cp);
+ if (!INTUSE(xdr_u_int) (xdrs, &u))
+ {
+ return FALSE;
+ }
+ *cp = u;
+ return TRUE;
+}
+
+/*
+ * XDR booleans
+ */
+bool_t
+xdr_bool (XDR *xdrs, bool_t *bp)
+{
+ long lb;
+
+ switch (xdrs->x_op)
+ {
+ case XDR_ENCODE:
+ lb = *bp ? XDR_TRUE : XDR_FALSE;
+ return XDR_PUTLONG (xdrs, &lb);
+
+ case XDR_DECODE:
+ if (!XDR_GETLONG (xdrs, &lb))
+ {
+ return FALSE;
+ }
+ *bp = (lb == XDR_FALSE) ? FALSE : TRUE;
+ return TRUE;
+
+ case XDR_FREE:
+ return TRUE;
+ }
+ return FALSE;
+}
+INTDEF(xdr_bool)
+
+/*
+ * XDR enumerations
+ */
+bool_t
+xdr_enum (XDR *xdrs, enum_t *ep)
+{
+ enum sizecheck
+ {
+ SIZEVAL
+ }; /* used to find the size of an enum */
+
+ /*
+ * enums are treated as ints
+ */
+ if (sizeof (enum sizecheck) == 4)
+ {
+#if INT_MAX < LONG_MAX
+ long l;
+
+ switch (xdrs->x_op)
+ {
+ case XDR_ENCODE:
+ l = *ep;
+ return XDR_PUTLONG (xdrs, &l);
+
+ case XDR_DECODE:
+ if (!XDR_GETLONG (xdrs, &l))
+ {
+ return FALSE;
+ }
+ *ep = l;
+ case XDR_FREE:
+ return TRUE;
+
+ }
+ return FALSE;
+#else
+ return INTUSE(xdr_long) (xdrs, (long *) ep);
+#endif
+ }
+ else if (sizeof (enum sizecheck) == sizeof (short))
+ {
+ return INTUSE(xdr_short) (xdrs, (short *) ep);
+ }
+ else
+ {
+ return FALSE;
+ }
+}
+INTDEF(xdr_enum)
+
+/*
+ * XDR opaque data
+ * Allows the specification of a fixed size sequence of opaque bytes.
+ * cp points to the opaque object and cnt gives the byte length.
+ */
+bool_t
+xdr_opaque (XDR *xdrs, caddr_t cp, u_int cnt)
+{
+ u_int rndup;
+ static char crud[BYTES_PER_XDR_UNIT];
+
+ /*
+ * if no data we are done
+ */
+ if (cnt == 0)
+ return TRUE;
+
+ /*
+ * round byte count to full xdr units
+ */
+ rndup = cnt % BYTES_PER_XDR_UNIT;
+ if (rndup > 0)
+ rndup = BYTES_PER_XDR_UNIT - rndup;
+
+ switch (xdrs->x_op)
+ {
+ case XDR_DECODE:
+ if (!XDR_GETBYTES (xdrs, cp, cnt))
+ {
+ return FALSE;
+ }
+ if (rndup == 0)
+ return TRUE;
+ return XDR_GETBYTES (xdrs, (caddr_t)crud, rndup);
+
+ case XDR_ENCODE:
+ if (!XDR_PUTBYTES (xdrs, cp, cnt))
+ {
+ return FALSE;
+ }
+ if (rndup == 0)
+ return TRUE;
+ return XDR_PUTBYTES (xdrs, xdr_zero, rndup);
+
+ case XDR_FREE:
+ return TRUE;
+ }
+ return FALSE;
+}
+INTDEF(xdr_opaque)
+
+/*
+ * XDR counted bytes
+ * *cpp is a pointer to the bytes, *sizep is the count.
+ * If *cpp is NULL maxsize bytes are allocated
+ */
+bool_t
+xdr_bytes (xdrs, cpp, sizep, maxsize)
+ XDR *xdrs;
+ char **cpp;
+ u_int *sizep;
+ u_int maxsize;
+{
+ char *sp = *cpp; /* sp is the actual string pointer */
+ u_int nodesize;
+
+ /*
+ * first deal with the length since xdr bytes are counted
+ */
+ if (!INTUSE(xdr_u_int) (xdrs, sizep))
+ {
+ return FALSE;
+ }
+ nodesize = *sizep;
+ if ((nodesize > maxsize) && (xdrs->x_op != XDR_FREE))
+ {
+ return FALSE;
+ }
+
+ /*
+ * now deal with the actual bytes
+ */
+ switch (xdrs->x_op)
+ {
+ case XDR_DECODE:
+ if (nodesize == 0)
+ {
+ return TRUE;
+ }
+ if (sp == NULL)
+ {
+ *cpp = sp = (char *) mem_alloc (nodesize);
+ }
+ if (sp == NULL)
+ {
+ fprintf (NULL, "%s", "xdr_bytes: out of memory\n");
+ return FALSE;
+ }
+ /* fall into ... */
+
+ case XDR_ENCODE:
+ return INTUSE(xdr_opaque) (xdrs, sp, nodesize);
+
+ case XDR_FREE:
+ if (sp != NULL)
+ {
+ mem_free (sp, nodesize);
+ *cpp = NULL;
+ }
+ return TRUE;
+ }
+ return FALSE;
+}
+INTDEF(xdr_bytes)
+
+/*
+ * Implemented here due to commonality of the object.
+ */
+bool_t
+xdr_netobj (xdrs, np)
+ XDR *xdrs;
+ struct netobj *np;
+{
+
+ return INTUSE(xdr_bytes) (xdrs, &np->n_bytes, &np->n_len, MAX_NETOBJ_SZ);
+}
+INTDEF(xdr_netobj)
+
+/*
+ * XDR a discriminated union
+ * Support routine for discriminated unions.
+ * You create an array of xdrdiscrim structures, terminated with
+ * an entry with a null procedure pointer. The routine gets
+ * the discriminant value and then searches the array of xdrdiscrims
+ * looking for that value. It calls the procedure given in the xdrdiscrim
+ * to handle the discriminant. If there is no specific routine a default
+ * routine may be called.
+ * If there is no specific or default routine an error is returned.
+ */
+bool_t
+xdr_union (xdrs, dscmp, unp, choices, dfault)
+ XDR *xdrs;
+ enum_t *dscmp; /* enum to decide which arm to work on */
+ char *unp; /* the union itself */
+ const struct xdr_discrim *choices; /* [value, xdr proc] for each arm */
+ xdrproc_t dfault; /* default xdr routine */
+{
+ enum_t dscm;
+
+ /*
+ * we deal with the discriminator; it's an enum
+ */
+ if (!INTUSE(xdr_enum) (xdrs, dscmp))
+ {
+ return FALSE;
+ }
+ dscm = *dscmp;
+
+ /*
+ * search choices for a value that matches the discriminator.
+ * if we find one, execute the xdr routine for that value.
+ */
+ for (; choices->proc != NULL_xdrproc_t; choices++)
+ {
+ if (choices->value == dscm)
+ return (*(choices->proc)) (xdrs, unp, LASTUNSIGNED);
+ }
+
+ /*
+ * no match - execute the default xdr routine if there is one
+ */
+ return ((dfault == NULL_xdrproc_t) ? FALSE :
+ (*dfault) (xdrs, unp, LASTUNSIGNED));
+}
+INTDEF(xdr_union)
+
+
+/*
+ * Non-portable xdr primitives.
+ * Care should be taken when moving these routines to new architectures.
+ */
+
+
+/*
+ * XDR null terminated ASCII strings
+ * xdr_string deals with "C strings" - arrays of bytes that are
+ * terminated by a NULL character. The parameter cpp references a
+ * pointer to storage; If the pointer is null, then the necessary
+ * storage is allocated. The last parameter is the max allowed length
+ * of the string as specified by a protocol.
+ */
+bool_t
+xdr_string (xdrs, cpp, maxsize)
+ XDR *xdrs;
+ char **cpp;
+ u_int maxsize;
+{
+ char *sp = *cpp; /* sp is the actual string pointer */
+ u_int size;
+ u_int nodesize;
+
+ /*
+ * first deal with the length since xdr strings are counted-strings
+ */
+ switch (xdrs->x_op)
+ {
+ case XDR_FREE:
+ if (sp == NULL)
+ {
+ return TRUE; /* already free */
+ }
+ /* fall through... */
+ case XDR_ENCODE:
+ if (sp == NULL)
+ return FALSE;
+ size = strlen (sp);
+ break;
+ case XDR_DECODE:
+ break;
+ }
+ if (!INTUSE(xdr_u_int) (xdrs, &size))
+ {
+ return FALSE;
+ }
+ if (size > maxsize)
+ {
+ return FALSE;
+ }
+ nodesize = size + 1;
+ if (nodesize == 0)
+ {
+ /* This means an overflow. It a bug in the caller which
+ provided a too large maxsize but nevertheless catch it
+ here. */
+ return FALSE;
+ }
+
+ /*
+ * now deal with the actual bytes
+ */
+ switch (xdrs->x_op)
+ {
+ case XDR_DECODE:
+ if (sp == NULL)
+ *cpp = sp = (char *) mem_alloc (nodesize);
+ if (sp == NULL)
+ {
+ fprintf (NULL, "%s", "xdr_string: out of memory\n");
+ return FALSE;
+ }
+ sp[size] = 0;
+ /* fall into ... */
+
+ case XDR_ENCODE:
+ return INTUSE(xdr_opaque) (xdrs, sp, size);
+
+ case XDR_FREE:
+ mem_free (sp, nodesize);
+ *cpp = NULL;
+ return TRUE;
+ }
+ return FALSE;
+}
+INTDEF(xdr_string)
+
+/*
+ * Wrapper for xdr_string that can be called directly from
+ * routines like clnt_call
+ */
+bool_t
+xdr_wrapstring (xdrs, cpp)
+ XDR *xdrs;
+ char **cpp;
+{
+ if (INTUSE(xdr_string) (xdrs, cpp, LASTUNSIGNED))
+ {
+ return TRUE;
+ }
+ return FALSE;
+}
--- /dev/null
+/*
+ * Sun RPC is a product of Sun Microsystems, Inc. and is provided for
+ * unrestricted use provided that this legend is included on all tape
+ * media and as a part of the software program in whole or part. Users
+ * may copy or modify Sun RPC without charge, but are not authorized
+ * to license or distribute it to anyone else except as part of a product or
+ * program developed by the user.
+ *
+ * SUN RPC IS PROVIDED AS IS WITH NO WARRANTIES OF ANY KIND INCLUDING THE
+ * WARRANTIES OF DESIGN, MERCHANTIBILITY AND FITNESS FOR A PARTICULAR
+ * PURPOSE, OR ARISING FROM A COURSE OF DEALING, USAGE OR TRADE PRACTICE.
+ *
+ * Sun RPC is provided with no support and without any obligation on the
+ * part of Sun Microsystems, Inc. to assist in its use, correction,
+ * modification or enhancement.
+ *
+ * SUN MICROSYSTEMS, INC. SHALL HAVE NO LIABILITY WITH RESPECT TO THE
+ * INFRINGEMENT OF COPYRIGHTS, TRADE SECRETS OR ANY PATENTS BY SUN RPC
+ * OR ANY PART THEREOF.
+ *
+ * In no event will Sun Microsystems, Inc. be liable for any lost revenue
+ * or profits or other special, indirect and consequential damages, even if
+ * Sun has been advised of the possibility of such damages.
+ *
+ * Sun Microsystems, Inc.
+ * 2550 Garcia Avenue
+ * Mountain View, California 94043
+ */
+
+/*
+ * xdr.h, External Data Representation Serialization Routines.
+ *
+ * Copyright (C) 1984, Sun Microsystems, Inc.
+ */
+
+#ifndef _RPC_XDR_H
+#define _RPC_XDR_H 1
+
+#include <features.h>
+#include <sys/types.h>
+#include "types.h"
+
+/* We need FILE. */
+#include <stdio.h>
+
+__BEGIN_DECLS
+
+/*
+ * XDR provides a conventional way for converting between C data
+ * types and an external bit-string representation. Library supplied
+ * routines provide for the conversion on built-in C data types. These
+ * routines and utility routines defined here are used to help implement
+ * a type encode/decode routine for each user-defined type.
+ *
+ * Each data type provides a single procedure which takes two arguments:
+ *
+ * bool_t
+ * xdrproc(xdrs, argresp)
+ * XDR *xdrs;
+ * <type> *argresp;
+ *
+ * xdrs is an instance of a XDR handle, to which or from which the data
+ * type is to be converted. argresp is a pointer to the structure to be
+ * converted. The XDR handle contains an operation field which indicates
+ * which of the operations (ENCODE, DECODE * or FREE) is to be performed.
+ *
+ * XDR_DECODE may allocate space if the pointer argresp is null. This
+ * data can be freed with the XDR_FREE operation.
+ *
+ * We write only one procedure per data type to make it easy
+ * to keep the encode and decode procedures for a data type consistent.
+ * In many cases the same code performs all operations on a user defined type,
+ * because all the hard work is done in the component type routines.
+ * decode as a series of calls on the nested data types.
+ */
+
+/*
+ * Xdr operations. XDR_ENCODE causes the type to be encoded into the
+ * stream. XDR_DECODE causes the type to be extracted from the stream.
+ * XDR_FREE can be used to release the space allocated by an XDR_DECODE
+ * request.
+ */
+enum xdr_op {
+ XDR_ENCODE = 0,
+ XDR_DECODE = 1,
+ XDR_FREE = 2
+};
+
+/*
+ * This is the number of bytes per unit of external data.
+ */
+#define BYTES_PER_XDR_UNIT (4)
+/*
+ * This only works if the above is a power of 2. But it's defined to be
+ * 4 by the appropriate RFCs. So it will work. And it's normally quicker
+ * than the old routine.
+ */
+#if 1
+#define RNDUP(x) (((x) + BYTES_PER_XDR_UNIT - 1) & ~(BYTES_PER_XDR_UNIT - 1))
+#else /* this is the old routine */
+#define RNDUP(x) ((((x) + BYTES_PER_XDR_UNIT - 1) / BYTES_PER_XDR_UNIT) \
+ * BYTES_PER_XDR_UNIT)
+#endif
+
+/*
+ * The XDR handle.
+ * Contains operation which is being applied to the stream,
+ * an operations vector for the particular implementation (e.g. see xdr_mem.c),
+ * and two private fields for the use of the particular implementation.
+ */
+typedef struct XDR XDR;
+struct XDR
+ {
+ enum xdr_op x_op; /* operation; fast additional param */
+ struct xdr_ops
+ {
+ bool_t (*x_getlong) (XDR *__xdrs, long *__lp);
+ /* get a long from underlying stream */
+ bool_t (*x_putlong) (XDR *__xdrs, __const long *__lp);
+ /* put a long to " */
+ bool_t (*x_getbytes) (XDR *__xdrs, caddr_t __addr, u_int __len);
+ /* get some bytes from " */
+ bool_t (*x_putbytes) (XDR *__xdrs, __const char *__addr, u_int __len);
+ /* put some bytes to " */
+ u_int (*x_getpostn) (__const XDR *__xdrs);
+ /* returns bytes off from beginning */
+ bool_t (*x_setpostn) (XDR *__xdrs, u_int __pos);
+ /* lets you reposition the stream */
+ int32_t *(*x_inline) (XDR *__xdrs, u_int __len);
+ /* buf quick ptr to buffered data */
+ void (*x_destroy) (XDR *__xdrs);
+ /* free privates of this xdr_stream */
+ bool_t (*x_getint32) (XDR *__xdrs, int32_t *__ip);
+ /* get a int from underlying stream */
+ bool_t (*x_putint32) (XDR *__xdrs, __const int32_t *__ip);
+ /* put a int to " */
+ }
+ *x_ops;
+ caddr_t x_public; /* users' data */
+ caddr_t x_private; /* pointer to private data */
+ caddr_t x_base; /* private used for position info */
+ u_int x_handy; /* extra private word */
+ };
+
+/*
+ * A xdrproc_t exists for each data type which is to be encoded or decoded.
+ *
+ * The second argument to the xdrproc_t is a pointer to an opaque pointer.
+ * The opaque pointer generally points to a structure of the data type
+ * to be decoded. If this pointer is 0, then the type routines should
+ * allocate dynamic storage of the appropriate size and return it.
+ * bool_t (*xdrproc_t)(XDR *, caddr_t *);
+ */
+typedef bool_t (*xdrproc_t) (XDR *, void *,...);
+
+
+/*
+ * Operations defined on a XDR handle
+ *
+ * XDR *xdrs;
+ * int32_t *int32p;
+ * long *longp;
+ * caddr_t addr;
+ * u_int len;
+ * u_int pos;
+ */
+#define XDR_GETINT32(xdrs, int32p) \
+ (*(xdrs)->x_ops->x_getint32)(xdrs, int32p)
+#define xdr_getint32(xdrs, int32p) \
+ (*(xdrs)->x_ops->x_getint32)(xdrs, int32p)
+
+#define XDR_PUTINT32(xdrs, int32p) \
+ (*(xdrs)->x_ops->x_putint32)(xdrs, int32p)
+#define xdr_putint32(xdrs, int32p) \
+ (*(xdrs)->x_ops->x_putint32)(xdrs, int32p)
+
+#define XDR_GETLONG(xdrs, longp) \
+ (*(xdrs)->x_ops->x_getlong)(xdrs, longp)
+#define xdr_getlong(xdrs, longp) \
+ (*(xdrs)->x_ops->x_getlong)(xdrs, longp)
+
+#define XDR_PUTLONG(xdrs, longp) \
+ (*(xdrs)->x_ops->x_putlong)(xdrs, longp)
+#define xdr_putlong(xdrs, longp) \
+ (*(xdrs)->x_ops->x_putlong)(xdrs, longp)
+
+#define XDR_GETBYTES(xdrs, addr, len) \
+ (*(xdrs)->x_ops->x_getbytes)(xdrs, addr, len)
+#define xdr_getbytes(xdrs, addr, len) \
+ (*(xdrs)->x_ops->x_getbytes)(xdrs, addr, len)
+
+#define XDR_PUTBYTES(xdrs, addr, len) \
+ (*(xdrs)->x_ops->x_putbytes)(xdrs, addr, len)
+#define xdr_putbytes(xdrs, addr, len) \
+ (*(xdrs)->x_ops->x_putbytes)(xdrs, addr, len)
+
+#define XDR_GETPOS(xdrs) \
+ (*(xdrs)->x_ops->x_getpostn)(xdrs)
+#define xdr_getpos(xdrs) \
+ (*(xdrs)->x_ops->x_getpostn)(xdrs)
+
+#define XDR_SETPOS(xdrs, pos) \
+ (*(xdrs)->x_ops->x_setpostn)(xdrs, pos)
+#define xdr_setpos(xdrs, pos) \
+ (*(xdrs)->x_ops->x_setpostn)(xdrs, pos)
+
+#define XDR_INLINE(xdrs, len) \
+ (*(xdrs)->x_ops->x_inline)(xdrs, len)
+#define xdr_inline(xdrs, len) \
+ (*(xdrs)->x_ops->x_inline)(xdrs, len)
+
+#define XDR_DESTROY(xdrs) \
+ do { \
+ if ((xdrs)->x_ops->x_destroy) \
+ (*(xdrs)->x_ops->x_destroy)(xdrs); \
+ } while (0)
+#define xdr_destroy(xdrs) \
+ do { \
+ if ((xdrs)->x_ops->x_destroy) \
+ (*(xdrs)->x_ops->x_destroy)(xdrs); \
+ } while (0)
+
+/*
+ * Support struct for discriminated unions.
+ * You create an array of xdrdiscrim structures, terminated with
+ * a entry with a null procedure pointer. The xdr_union routine gets
+ * the discriminant value and then searches the array of structures
+ * for a matching value. If a match is found the associated xdr routine
+ * is called to handle that part of the union. If there is
+ * no match, then a default routine may be called.
+ * If there is no match and no default routine it is an error.
+ */
+#define NULL_xdrproc_t ((xdrproc_t)0)
+struct xdr_discrim
+{
+ int value;
+ xdrproc_t proc;
+};
+
+/*
+ * Inline routines for fast encode/decode of primitive data types.
+ * Caveat emptor: these use single memory cycles to get the
+ * data from the underlying buffer, and will fail to operate
+ * properly if the data is not aligned. The standard way to use these
+ * is to say:
+ * if ((buf = XDR_INLINE(xdrs, count)) == NULL)
+ * return (FALSE);
+ * <<< macro calls >>>
+ * where ``count'' is the number of bytes of data occupied
+ * by the primitive data types.
+ *
+ * N.B. and frozen for all time: each data type here uses 4 bytes
+ * of external representation.
+ */
+
+#define IXDR_GET_INT32(buf) ((int32_t)ntohl((uint32_t)*(buf)++))
+#define IXDR_PUT_INT32(buf, v) (*(buf)++ = (int32_t)htonl((uint32_t)(v)))
+#define IXDR_GET_U_INT32(buf) ((uint32_t)IXDR_GET_INT32(buf))
+#define IXDR_PUT_U_INT32(buf, v) IXDR_PUT_INT32(buf, (int32_t)(v))
+
+/* WARNING: The IXDR_*_LONG defines are removed by Sun for new platforms
+ * and shouldn't be used any longer. Code which use this defines or longs
+ * in the RPC code will not work on 64bit Solaris platforms !
+ */
+#define IXDR_GET_LONG(buf) ((long)IXDR_GET_U_INT32(buf))
+#define IXDR_PUT_LONG(buf, v) ((long)IXDR_PUT_INT32(buf, (long)(v)))
+#define IXDR_GET_U_LONG(buf) ((u_long)IXDR_GET_LONG(buf))
+#define IXDR_PUT_U_LONG(buf, v) IXDR_PUT_LONG(buf, (long)(v))
+
+
+#define IXDR_GET_BOOL(buf) ((bool_t)IXDR_GET_LONG(buf))
+#define IXDR_GET_ENUM(buf, t) ((t)IXDR_GET_LONG(buf))
+#define IXDR_GET_SHORT(buf) ((short)IXDR_GET_LONG(buf))
+#define IXDR_GET_U_SHORT(buf) ((u_short)IXDR_GET_LONG(buf))
+
+#define IXDR_PUT_BOOL(buf, v) IXDR_PUT_LONG(buf, (long)(v))
+#define IXDR_PUT_ENUM(buf, v) IXDR_PUT_LONG(buf, (long)(v))
+#define IXDR_PUT_SHORT(buf, v) IXDR_PUT_LONG(buf, (long)(v))
+#define IXDR_PUT_U_SHORT(buf, v) IXDR_PUT_LONG(buf, (long)(v))
+
+/*
+ * These are the "generic" xdr routines.
+ * None of these can have const applied because it's not possible to
+ * know whether the call is a read or a write to the passed parameter
+ * also, the XDR structure is always updated by some of these calls.
+ */
+extern bool_t xdr_void (void) __THROW;
+extern bool_t xdr_short (XDR *__xdrs, short *__sp) __THROW;
+extern bool_t xdr_u_short (XDR *__xdrs, u_short *__usp) __THROW;
+extern bool_t xdr_int (XDR *__xdrs, int *__ip) __THROW;
+extern bool_t xdr_u_int (XDR *__xdrs, u_int *__up) __THROW;
+extern bool_t xdr_long (XDR *__xdrs, long *__lp) __THROW;
+extern bool_t xdr_u_long (XDR *__xdrs, u_long *__ulp) __THROW;
+extern bool_t xdr_hyper (XDR *__xdrs, quad_t *__llp) __THROW;
+extern bool_t xdr_u_hyper (XDR *__xdrs, u_quad_t *__ullp) __THROW;
+extern bool_t xdr_longlong_t (XDR *__xdrs, quad_t *__llp) __THROW;
+extern bool_t xdr_u_longlong_t (XDR *__xdrs, u_quad_t *__ullp) __THROW;
+extern bool_t xdr_int8_t (XDR *__xdrs, int8_t *__ip) __THROW;
+extern bool_t xdr_uint8_t (XDR *__xdrs, uint8_t *__up) __THROW;
+extern bool_t xdr_int16_t (XDR *__xdrs, int16_t *__ip) __THROW;
+extern bool_t xdr_uint16_t (XDR *__xdrs, uint16_t *__up) __THROW;
+extern bool_t xdr_int32_t (XDR *__xdrs, int32_t *__ip) __THROW;
+extern bool_t xdr_uint32_t (XDR *__xdrs, uint32_t *__up) __THROW;
+extern bool_t xdr_int64_t (XDR *__xdrs, int64_t *__ip) __THROW;
+extern bool_t xdr_uint64_t (XDR *__xdrs, uint64_t *__up) __THROW;
+extern bool_t xdr_quad_t (XDR *__xdrs, quad_t *__ip) __THROW;
+extern bool_t xdr_u_quad_t (XDR *__xdrs, u_quad_t *__up) __THROW;
+extern bool_t xdr_bool (XDR *__xdrs, bool_t *__bp) __THROW;
+extern bool_t xdr_enum (XDR *__xdrs, enum_t *__ep) __THROW;
+extern bool_t xdr_array (XDR * _xdrs, caddr_t *__addrp, u_int *__sizep,
+ u_int __maxsize, u_int __elsize, xdrproc_t __elproc)
+ __THROW;
+extern bool_t xdr_bytes (XDR *__xdrs, char **__cpp, u_int *__sizep,
+ u_int __maxsize) __THROW;
+extern bool_t xdr_opaque (XDR *__xdrs, caddr_t __cp, u_int __cnt) __THROW;
+extern bool_t xdr_string (XDR *__xdrs, char **__cpp, u_int __maxsize) __THROW;
+extern bool_t xdr_union (XDR *__xdrs, enum_t *__dscmp, char *__unp,
+ __const struct xdr_discrim *__choices,
+ xdrproc_t dfault) __THROW;
+extern bool_t xdr_char (XDR *__xdrs, char *__cp) __THROW;
+extern bool_t xdr_u_char (XDR *__xdrs, u_char *__cp) __THROW;
+extern bool_t xdr_vector (XDR *__xdrs, char *__basep, u_int __nelem,
+ u_int __elemsize, xdrproc_t __xdr_elem) __THROW;
+extern bool_t xdr_float (XDR *__xdrs, float *__fp) __THROW;
+extern bool_t xdr_double (XDR *__xdrs, double *__dp) __THROW;
+extern bool_t xdr_reference (XDR *__xdrs, caddr_t *__xpp, u_int __size,
+ xdrproc_t __proc) __THROW;
+extern bool_t xdr_pointer (XDR *__xdrs, char **__objpp,
+ u_int __obj_size, xdrproc_t __xdr_obj) __THROW;
+extern bool_t xdr_wrapstring (XDR *__xdrs, char **__cpp) __THROW;
+extern u_long xdr_sizeof (xdrproc_t, void *) __THROW;
+
+/*
+ * Common opaque bytes objects used by many rpc protocols;
+ * declared here due to commonality.
+ */
+#define MAX_NETOBJ_SZ 1024
+struct netobj
+{
+ u_int n_len;
+ char *n_bytes;
+};
+typedef struct netobj netobj;
+extern bool_t xdr_netobj (XDR *__xdrs, struct netobj *__np) __THROW;
+
+/*
+ * These are the public routines for the various implementations of
+ * xdr streams.
+ */
+
+/* XDR using memory buffers */
+extern void xdrmem_create (XDR *__xdrs, __const caddr_t __addr,
+ u_int __size, enum xdr_op __xop) __THROW;
+
+/* XDR using stdio library */
+extern void xdrstdio_create (XDR *__xdrs, FILE *__file, enum xdr_op __xop)
+ __THROW;
+
+/* XDR pseudo records for tcp */
+extern void xdrrec_create (XDR *__xdrs, u_int __sendsize,
+ u_int __recvsize, caddr_t __tcp_handle,
+ int (*__readit) (char *, char *, int),
+ int (*__writeit) (char *, char *, int)) __THROW;
+
+/* make end of xdr record */
+extern bool_t xdrrec_endofrecord (XDR *__xdrs, bool_t __sendnow) __THROW;
+
+/* move to beginning of next record */
+extern bool_t xdrrec_skiprecord (XDR *__xdrs) __THROW;
+
+/* true if no more input */
+extern bool_t xdrrec_eof (XDR *__xdrs) __THROW;
+
+/* free memory buffers for xdr */
+extern void xdr_free (xdrproc_t __proc, char *__objp) __THROW;
+
+__END_DECLS
+
+#endif /* rpc/xdr.h */
--- /dev/null
+# define INTUSE(name) name
+# define INTDEF(name)
+/* @(#)xdr_array.c 2.1 88/07/29 4.0 RPCSRC */
+/*
+ * Sun RPC is a product of Sun Microsystems, Inc. and is provided for
+ * unrestricted use provided that this legend is included on all tape
+ * media and as a part of the software program in whole or part. Users
+ * may copy or modify Sun RPC without charge, but are not authorized
+ * to license or distribute it to anyone else except as part of a product or
+ * program developed by the user.
+ *
+ * SUN RPC IS PROVIDED AS IS WITH NO WARRANTIES OF ANY KIND INCLUDING THE
+ * WARRANTIES OF DESIGN, MERCHANTIBILITY AND FITNESS FOR A PARTICULAR
+ * PURPOSE, OR ARISING FROM A COURSE OF DEALING, USAGE OR TRADE PRACTICE.
+ *
+ * Sun RPC is provided with no support and without any obligation on the
+ * part of Sun Microsystems, Inc. to assist in its use, correction,
+ * modification or enhancement.
+ *
+ * SUN MICROSYSTEMS, INC. SHALL HAVE NO LIABILITY WITH RESPECT TO THE
+ * INFRINGEMENT OF COPYRIGHTS, TRADE SECRETS OR ANY PATENTS BY SUN RPC
+ * OR ANY PART THEREOF.
+ *
+ * In no event will Sun Microsystems, Inc. be liable for any lost revenue
+ * or profits or other special, indirect and consequential damages, even if
+ * Sun has been advised of the possibility of such damages.
+ *
+ * Sun Microsystems, Inc.
+ * 2550 Garcia Avenue
+ * Mountain View, California 94043
+ */
+#if !defined(lint) && defined(SCCSIDS)
+static char sccsid[] = "@(#)xdr_array.c 1.10 87/08/11 Copyr 1984 Sun Micro";
+#endif
+
+/*
+ * xdr_array.c, Generic XDR routines implementation.
+ *
+ * Copyright (C) 1984, Sun Microsystems, Inc.
+ *
+ * These are the "non-trivial" xdr primitives used to serialize and de-serialize
+ * arrays. See xdr.h for more info on the interface to xdr.
+ */
+
+#include <stdio.h>
+#include <string.h>
+#include "types.h"
+#include "xdr.h"
+#include <libintl.h>
+#include <limits.h>
+
+#ifdef USE_IN_LIBIO
+# include <wchar.h>
+#endif
+
+#define LASTUNSIGNED ((u_int)0-1)
+
+
+/*
+ * XDR an array of arbitrary elements
+ * *addrp is a pointer to the array, *sizep is the number of elements.
+ * If addrp is NULL (*sizep * elsize) bytes are allocated.
+ * elsize is the size (in bytes) of each element, and elproc is the
+ * xdr procedure to call to handle each element of the array.
+ */
+bool_t
+xdr_array (xdrs, addrp, sizep, maxsize, elsize, elproc)
+ XDR *xdrs;
+ caddr_t *addrp; /* array pointer */
+ u_int *sizep; /* number of elements */
+ u_int maxsize; /* max numberof elements */
+ u_int elsize; /* size in bytes of each element */
+ xdrproc_t elproc; /* xdr routine to handle each element */
+{
+ u_int i;
+ caddr_t target = *addrp;
+ u_int c; /* the actual element count */
+ bool_t stat = TRUE;
+ u_int nodesize;
+
+ /* like strings, arrays are really counted arrays */
+ if (!INTUSE(xdr_u_int) (xdrs, sizep))
+ {
+ return FALSE;
+ }
+ c = *sizep;
+ /*
+ * XXX: Let the overflow possibly happen with XDR_FREE because mem_free()
+ * doesn't actually use its second argument anyway.
+ */
+ if ((c > maxsize || c > UINT_MAX / elsize) && (xdrs->x_op != XDR_FREE))
+ {
+ return FALSE;
+ }
+ nodesize = c * elsize;
+
+ /*
+ * if we are deserializing, we may need to allocate an array.
+ * We also save time by checking for a null array if we are freeing.
+ */
+ if (target == NULL)
+ switch (xdrs->x_op)
+ {
+ case XDR_DECODE:
+ if (c == 0)
+ return TRUE;
+ *addrp = target = mem_alloc (nodesize);
+ if (target == NULL)
+ {
+ fprintf (stderr, "%s", "xdr_array: out of memory\n");
+ return FALSE;
+ }
+ __bzero (target, nodesize);
+ break;
+
+ case XDR_FREE:
+ return TRUE;
+ default:
+ break;
+ }
+
+ /*
+ * now we xdr each element of array
+ */
+ for (i = 0; (i < c) && stat; i++)
+ {
+ stat = (*elproc) (xdrs, target, LASTUNSIGNED);
+ target += elsize;
+ }
+
+ /*
+ * the array may need freeing
+ */
+ if (xdrs->x_op == XDR_FREE)
+ {
+ mem_free (*addrp, nodesize);
+ *addrp = NULL;
+ }
+ return stat;
+}
+INTDEF(xdr_array)
+
+/*
+ * xdr_vector():
+ *
+ * XDR a fixed length array. Unlike variable-length arrays,
+ * the storage of fixed length arrays is static and unfreeable.
+ * > basep: base of the array
+ * > size: size of the array
+ * > elemsize: size of each element
+ * > xdr_elem: routine to XDR each element
+ */
+bool_t
+xdr_vector (xdrs, basep, nelem, elemsize, xdr_elem)
+ XDR *xdrs;
+ char *basep;
+ u_int nelem;
+ u_int elemsize;
+ xdrproc_t xdr_elem;
+{
+ u_int i;
+ char *elptr;
+
+ elptr = basep;
+ for (i = 0; i < nelem; i++)
+ {
+ if (!(*xdr_elem) (xdrs, elptr, LASTUNSIGNED))
+ {
+ return FALSE;
+ }
+ elptr += elemsize;
+ }
+ return TRUE;
+}
--- /dev/null
+/* @(#)xdr_float.c 2.1 88/07/29 4.0 RPCSRC */
+/*
+ * Sun RPC is a product of Sun Microsystems, Inc. and is provided for
+ * unrestricted use provided that this legend is included on all tape
+ * media and as a part of the software program in whole or part. Users
+ * may copy or modify Sun RPC without charge, but are not authorized
+ * to license or distribute it to anyone else except as part of a product or
+ * program developed by the user.
+ *
+ * SUN RPC IS PROVIDED AS IS WITH NO WARRANTIES OF ANY KIND INCLUDING THE
+ * WARRANTIES OF DESIGN, MERCHANTIBILITY AND FITNESS FOR A PARTICULAR
+ * PURPOSE, OR ARISING FROM A COURSE OF DEALING, USAGE OR TRADE PRACTICE.
+ *
+ * Sun RPC is provided with no support and without any obligation on the
+ * part of Sun Microsystems, Inc. to assist in its use, correction,
+ * modification or enhancement.
+ *
+ * SUN MICROSYSTEMS, INC. SHALL HAVE NO LIABILITY WITH RESPECT TO THE
+ * INFRINGEMENT OF COPYRIGHTS, TRADE SECRETS OR ANY PATENTS BY SUN RPC
+ * OR ANY PART THEREOF.
+ *
+ * In no event will Sun Microsystems, Inc. be liable for any lost revenue
+ * or profits or other special, indirect and consequential damages, even if
+ * Sun has been advised of the possibility of such damages.
+ *
+ * Sun Microsystems, Inc.
+ * 2550 Garcia Avenue
+ * Mountain View, California 94043
+ */
+#if !defined(lint) && defined(SCCSIDS)
+static char sccsid[] = "@(#)xdr_float.c 1.12 87/08/11 Copyr 1984 Sun Micro";
+#endif
+
+/*
+ * xdr_float.c, Generic XDR routines implementation.
+ *
+ * Copyright (C) 1984, Sun Microsystems, Inc.
+ *
+ * These are the "floating point" xdr routines used to (de)serialize
+ * most common data items. See xdr.h for more info on the interface to
+ * xdr.
+ */
+
+#include <stdio.h>
+#include <endian.h>
+
+#include "types.h"
+#include "xdr.h"
+
+/*
+ * NB: Not portable.
+ * This routine works on Suns (Sky / 68000's) and Vaxen.
+ */
+
+#define LSW (__FLOAT_WORD_ORDER == __BIG_ENDIAN)
+
+#ifdef vax
+
+/* What IEEE single precision floating point looks like on a Vax */
+struct ieee_single {
+ unsigned int mantissa: 23;
+ unsigned int exp : 8;
+ unsigned int sign : 1;
+};
+
+/* Vax single precision floating point */
+struct vax_single {
+ unsigned int mantissa1 : 7;
+ unsigned int exp : 8;
+ unsigned int sign : 1;
+ unsigned int mantissa2 : 16;
+};
+
+#define VAX_SNG_BIAS 0x81
+#define IEEE_SNG_BIAS 0x7f
+
+static struct sgl_limits {
+ struct vax_single s;
+ struct ieee_single ieee;
+} sgl_limits[2] = {
+ {{ 0x7f, 0xff, 0x0, 0xffff }, /* Max Vax */
+ { 0x0, 0xff, 0x0 }}, /* Max IEEE */
+ {{ 0x0, 0x0, 0x0, 0x0 }, /* Min Vax */
+ { 0x0, 0x0, 0x0 }} /* Min IEEE */
+};
+#endif /* vax */
+
+bool_t
+xdr_float(xdrs, fp)
+ XDR *xdrs;
+ float *fp;
+{
+#ifdef vax
+ struct ieee_single is;
+ struct vax_single vs, *vsp;
+ struct sgl_limits *lim;
+ int i;
+#endif
+ switch (xdrs->x_op) {
+
+ case XDR_ENCODE:
+#ifdef vax
+ vs = *((struct vax_single *)fp);
+ for (i = 0, lim = sgl_limits;
+ i < sizeof(sgl_limits)/sizeof(struct sgl_limits);
+ i++, lim++) {
+ if ((vs.mantissa2 == lim->s.mantissa2) &&
+ (vs.exp == lim->s.exp) &&
+ (vs.mantissa1 == lim->s.mantissa1)) {
+ is = lim->ieee;
+ goto shipit;
+ }
+ }
+ is.exp = vs.exp - VAX_SNG_BIAS + IEEE_SNG_BIAS;
+ is.mantissa = (vs.mantissa1 << 16) | vs.mantissa2;
+ shipit:
+ is.sign = vs.sign;
+ return (XDR_PUTLONG(xdrs, (long *)&is));
+#else
+ if (sizeof(float) == sizeof(long))
+ return (XDR_PUTLONG(xdrs, (long *)fp));
+ else if (sizeof(float) == sizeof(int)) {
+ long tmp = *(int *)fp;
+ return (XDR_PUTLONG(xdrs, &tmp));
+ }
+ break;
+#endif
+
+ case XDR_DECODE:
+#ifdef vax
+ vsp = (struct vax_single *)fp;
+ if (!XDR_GETLONG(xdrs, (long *)&is))
+ return (FALSE);
+ for (i = 0, lim = sgl_limits;
+ i < sizeof(sgl_limits)/sizeof(struct sgl_limits);
+ i++, lim++) {
+ if ((is.exp == lim->ieee.exp) &&
+ (is.mantissa == lim->ieee.mantissa)) {
+ *vsp = lim->s;
+ goto doneit;
+ }
+ }
+ vsp->exp = is.exp - IEEE_SNG_BIAS + VAX_SNG_BIAS;
+ vsp->mantissa2 = is.mantissa;
+ vsp->mantissa1 = (is.mantissa >> 16);
+ doneit:
+ vsp->sign = is.sign;
+ return (TRUE);
+#else
+ if (sizeof(float) == sizeof(long))
+ return (XDR_GETLONG(xdrs, (long *)fp));
+ else if (sizeof(float) == sizeof(int)) {
+ long tmp;
+ if (XDR_GETLONG(xdrs, &tmp)) {
+ *(int *)fp = tmp;
+ return (TRUE);
+ }
+ }
+ break;
+#endif
+
+ case XDR_FREE:
+ return (TRUE);
+ }
+ return (FALSE);
+}
+
+/*
+ * This routine works on Suns (Sky / 68000's) and Vaxen.
+ */
+
+#ifdef vax
+/* What IEEE double precision floating point looks like on a Vax */
+struct ieee_double {
+ unsigned int mantissa1 : 20;
+ unsigned int exp : 11;
+ unsigned int sign : 1;
+ unsigned int mantissa2 : 32;
+};
+
+/* Vax double precision floating point */
+struct vax_double {
+ unsigned int mantissa1 : 7;
+ unsigned int exp : 8;
+ unsigned int sign : 1;
+ unsigned int mantissa2 : 16;
+ unsigned int mantissa3 : 16;
+ unsigned int mantissa4 : 16;
+};
+
+#define VAX_DBL_BIAS 0x81
+#define IEEE_DBL_BIAS 0x3ff
+#define MASK(nbits) ((1 << nbits) - 1)
+
+static struct dbl_limits {
+ struct vax_double d;
+ struct ieee_double ieee;
+} dbl_limits[2] = {
+ {{ 0x7f, 0xff, 0x0, 0xffff, 0xffff, 0xffff }, /* Max Vax */
+ { 0x0, 0x7ff, 0x0, 0x0 }}, /* Max IEEE */
+ {{ 0x0, 0x0, 0x0, 0x0, 0x0, 0x0}, /* Min Vax */
+ { 0x0, 0x0, 0x0, 0x0 }} /* Min IEEE */
+};
+
+#endif /* vax */
+
+
+bool_t
+xdr_double(xdrs, dp)
+ XDR *xdrs;
+ double *dp;
+{
+#ifdef vax
+ struct ieee_double id;
+ struct vax_double vd;
+ register struct dbl_limits *lim;
+ int i;
+#endif
+
+ switch (xdrs->x_op) {
+
+ case XDR_ENCODE:
+#ifdef vax
+ vd = *((struct vax_double *)dp);
+ for (i = 0, lim = dbl_limits;
+ i < sizeof(dbl_limits)/sizeof(struct dbl_limits);
+ i++, lim++) {
+ if ((vd.mantissa4 == lim->d.mantissa4) &&
+ (vd.mantissa3 == lim->d.mantissa3) &&
+ (vd.mantissa2 == lim->d.mantissa2) &&
+ (vd.mantissa1 == lim->d.mantissa1) &&
+ (vd.exp == lim->d.exp)) {
+ id = lim->ieee;
+ goto shipit;
+ }
+ }
+ id.exp = vd.exp - VAX_DBL_BIAS + IEEE_DBL_BIAS;
+ id.mantissa1 = (vd.mantissa1 << 13) | (vd.mantissa2 >> 3);
+ id.mantissa2 = ((vd.mantissa2 & MASK(3)) << 29) |
+ (vd.mantissa3 << 13) |
+ ((vd.mantissa4 >> 3) & MASK(13));
+ shipit:
+ id.sign = vd.sign;
+ dp = (double *)&id;
+#endif
+ if (2*sizeof(long) == sizeof(double)) {
+ long *lp = (long *)dp;
+ return (XDR_PUTLONG(xdrs, lp+!LSW) &&
+ XDR_PUTLONG(xdrs, lp+LSW));
+ } else if (2*sizeof(int) == sizeof(double)) {
+ int *ip = (int *)dp;
+ long tmp[2];
+ tmp[0] = ip[!LSW];
+ tmp[1] = ip[LSW];
+ return (XDR_PUTLONG(xdrs, tmp) &&
+ XDR_PUTLONG(xdrs, tmp+1));
+ }
+ break;
+
+ case XDR_DECODE:
+#ifdef vax
+ lp = (long *)&id;
+ if (!XDR_GETLONG(xdrs, lp++) || !XDR_GETLONG(xdrs, lp))
+ return (FALSE);
+ for (i = 0, lim = dbl_limits;
+ i < sizeof(dbl_limits)/sizeof(struct dbl_limits);
+ i++, lim++) {
+ if ((id.mantissa2 == lim->ieee.mantissa2) &&
+ (id.mantissa1 == lim->ieee.mantissa1) &&
+ (id.exp == lim->ieee.exp)) {
+ vd = lim->d;
+ goto doneit;
+ }
+ }
+ vd.exp = id.exp - IEEE_DBL_BIAS + VAX_DBL_BIAS;
+ vd.mantissa1 = (id.mantissa1 >> 13);
+ vd.mantissa2 = ((id.mantissa1 & MASK(13)) << 3) |
+ (id.mantissa2 >> 29);
+ vd.mantissa3 = (id.mantissa2 >> 13);
+ vd.mantissa4 = (id.mantissa2 << 3);
+ doneit:
+ vd.sign = id.sign;
+ *dp = *((double *)&vd);
+ return (TRUE);
+#else
+ if (2*sizeof(long) == sizeof(double)) {
+ long *lp = (long *)dp;
+ return (XDR_GETLONG(xdrs, lp+!LSW) &&
+ XDR_GETLONG(xdrs, lp+LSW));
+ } else if (2*sizeof(int) == sizeof(double)) {
+ int *ip = (int *)dp;
+ long tmp[2];
+ if (XDR_GETLONG(xdrs, tmp+!LSW) &&
+ XDR_GETLONG(xdrs, tmp+LSW)) {
+ ip[0] = tmp[0];
+ ip[1] = tmp[1];
+ return (TRUE);
+ }
+ }
+ break;
+#endif
+
+ case XDR_FREE:
+ return (TRUE);
+ }
+ return (FALSE);
+}
--- /dev/null
+/*
+ * Sun RPC is a product of Sun Microsystems, Inc. and is provided for
+ * unrestricted use provided that this legend is included on all tape
+ * media and as a part of the software program in whole or part. Users
+ * may copy or modify Sun RPC without charge, but are not authorized
+ * to license or distribute it to anyone else except as part of a product or
+ * program developed by the user.
+ *
+ * SUN RPC IS PROVIDED AS IS WITH NO WARRANTIES OF ANY KIND INCLUDING THE
+ * WARRANTIES OF DESIGN, MERCHANTIBILITY AND FITNESS FOR A PARTICULAR
+ * PURPOSE, OR ARISING FROM A COURSE OF DEALING, USAGE OR TRADE PRACTICE.
+ *
+ * Sun RPC is provided with no support and without any obligation on the
+ * part of Sun Microsystems, Inc. to assist in its use, correction,
+ * modification or enhancement.
+ *
+ * SUN MICROSYSTEMS, INC. SHALL HAVE NO LIABILITY WITH RESPECT TO THE
+ * INFRINGEMENT OF COPYRIGHTS, TRADE SECRETS OR ANY PATENTS BY SUN RPC
+ * OR ANY PART THEREOF.
+ *
+ * In no event will Sun Microsystems, Inc. be liable for any lost revenue
+ * or profits or other special, indirect and consequential damages, even if
+ * Sun has been advised of the possibility of such damages.
+ *
+ * Sun Microsystems, Inc.
+ * 2550 Garcia Avenue
+ * Mountain View, California 94043
+ */
+
+/*
+ * xdr_stdio.c, XDR implementation on standard i/o file.
+ *
+ * Copyright (C) 1984, Sun Microsystems, Inc.
+ *
+ * This set of routines implements a XDR on a stdio stream.
+ * XDR_ENCODE serializes onto the stream, XDR_DECODE de-serializes
+ * from the stream.
+ */
+
+#include "types.h"
+#include <stdio.h>
+#include "xdr.h"
+
+#ifdef USE_IN_LIBIO
+# include <libio/iolibio.h>
+# define fflush(s) INTUSE(_IO_fflush) (s)
+# define fread(p, m, n, s) INTUSE(_IO_fread) (p, m, n, s)
+# define ftell(s) INTUSE(_IO_ftell) (s)
+# define fwrite(p, m, n, s) INTUSE(_IO_fwrite) (p, m, n, s)
+#endif
+
+static bool_t xdrstdio_getlong (XDR *, long *);
+static bool_t xdrstdio_putlong (XDR *, const long *);
+static bool_t xdrstdio_getbytes (XDR *, caddr_t, u_int);
+static bool_t xdrstdio_putbytes (XDR *, const char *, u_int);
+static u_int xdrstdio_getpos (const XDR *);
+static bool_t xdrstdio_setpos (XDR *, u_int);
+static int32_t *xdrstdio_inline (XDR *, u_int);
+static void xdrstdio_destroy (XDR *);
+static bool_t xdrstdio_getint32 (XDR *, int32_t *);
+static bool_t xdrstdio_putint32 (XDR *, const int32_t *);
+
+/*
+ * Ops vector for stdio type XDR
+ */
+static const struct xdr_ops xdrstdio_ops =
+{
+ xdrstdio_getlong, /* deserialize a long int */
+ xdrstdio_putlong, /* serialize a long int */
+ xdrstdio_getbytes, /* deserialize counted bytes */
+ xdrstdio_putbytes, /* serialize counted bytes */
+ xdrstdio_getpos, /* get offset in the stream */
+ xdrstdio_setpos, /* set offset in the stream */
+ xdrstdio_inline, /* prime stream for inline macros */
+ xdrstdio_destroy, /* destroy stream */
+ xdrstdio_getint32, /* deserialize a int */
+ xdrstdio_putint32 /* serialize a int */
+};
+
+/*
+ * Initialize a stdio xdr stream.
+ * Sets the xdr stream handle xdrs for use on the stream file.
+ * Operation flag is set to op.
+ */
+void
+xdrstdio_create (XDR *xdrs, FILE *file, enum xdr_op op)
+{
+ xdrs->x_op = op;
+ /* We have to add the const since the `struct xdr_ops' in `struct XDR'
+ is not `const'. */
+ xdrs->x_ops = (struct xdr_ops *) &xdrstdio_ops;
+ xdrs->x_private = (caddr_t) file;
+ xdrs->x_handy = 0;
+ xdrs->x_base = 0;
+}
+
+/*
+ * Destroy a stdio xdr stream.
+ * Cleans up the xdr stream handle xdrs previously set up by xdrstdio_create.
+ */
+static void
+xdrstdio_destroy (XDR *xdrs)
+{
+ (void) fflush ((FILE *) xdrs->x_private);
+ /* xx should we close the file ?? */
+};
+
+static bool_t
+xdrstdio_getlong (XDR *xdrs, long *lp)
+{
+ u_int32_t mycopy;
+
+ if (fread ((caddr_t) &mycopy, 4, 1, (FILE *) xdrs->x_private) != 1)
+ return FALSE;
+ *lp = (long) ntohl (mycopy);
+ return TRUE;
+}
+
+static bool_t
+xdrstdio_putlong (XDR *xdrs, const long *lp)
+{
+ int32_t mycopy = htonl ((u_int32_t) *lp);
+
+ if (fwrite ((caddr_t) &mycopy, 4, 1, (FILE *) xdrs->x_private) != 1)
+ return FALSE;
+ return TRUE;
+}
+
+static bool_t
+xdrstdio_getbytes (XDR *xdrs, const caddr_t addr, u_int len)
+{
+ if ((len != 0) && (fread (addr, (int) len, 1,
+ (FILE *) xdrs->x_private) != 1))
+ return FALSE;
+ return TRUE;
+}
+
+static bool_t
+xdrstdio_putbytes (XDR *xdrs, const char *addr, u_int len)
+{
+ if ((len != 0) && (fwrite (addr, (int) len, 1,
+ (FILE *) xdrs->x_private) != 1))
+ return FALSE;
+ return TRUE;
+}
+
+static u_int
+xdrstdio_getpos (const XDR *xdrs)
+{
+ return (u_int) ftell ((FILE *) xdrs->x_private);
+}
+
+static bool_t
+xdrstdio_setpos (XDR *xdrs, u_int pos)
+{
+ return fseek ((FILE *) xdrs->x_private, (long) pos, 0) < 0 ? FALSE : TRUE;
+}
+
+static int32_t *
+xdrstdio_inline (XDR *xdrs, u_int len)
+{
+ /*
+ * Must do some work to implement this: must insure
+ * enough data in the underlying stdio buffer,
+ * that the buffer is aligned so that we can indirect through a
+ * long *, and stuff this pointer in xdrs->x_buf. Doing
+ * a fread or fwrite to a scratch buffer would defeat
+ * most of the gains to be had here and require storage
+ * management on this buffer, so we don't do this.
+ */
+ return NULL;
+}
+
+static bool_t
+xdrstdio_getint32 (XDR *xdrs, int32_t *ip)
+{
+ int32_t mycopy;
+
+ if (fread ((caddr_t) &mycopy, 4, 1, (FILE *) xdrs->x_private) != 1)
+ return FALSE;
+ *ip = ntohl (mycopy);
+ return TRUE;
+}
+
+static bool_t
+xdrstdio_putint32 (XDR *xdrs, const int32_t *ip)
+{
+ int32_t mycopy = htonl (*ip);
+
+ ip = &mycopy;
+ if (fwrite ((caddr_t) ip, 4, 1, (FILE *) xdrs->x_private) != 1)
+ return FALSE;
+ return TRUE;
+}
+
+/* libc_hidden_def (xdrstdio_create) */
--- /dev/null
+/*_________________________________________________________________
+ |
+ | xdrf.h - include file for C routines that want to use the
+ | functions below.
+*/
+
+int xdropen(XDR *xdrs, const char *filename, const char *type);
+int xdrclose(XDR *xdrs) ;
+int xdr3dfcoord(XDR *xdrs, float *fp, int *size, float *precision) ;
+
integer iscode,indpdb,outpdb,outmol2,icomparfunc,pdbint,
- & ensembles
+ & ensembles,constr_dist
logical refstr,pdbref,punch_dist,print_rms,caonly,verbose,
& merge_helices,bxfile,cxfile,histfile,entfile,zscfile,
& rmsrgymap,with_dihed_constr,check_conf,histout
common /cntrl/ iscode,indpdb,refstr,pdbref,outpdb,outmol2,
& punch_dist,print_rms,caonly,verbose,icomparfunc,pdbint,
& merge_helices,bxfile,cxfile,histfile,entfile,zscfile,rmsrgymap,
- & ensembles,with_dihed_constr,check_conf,histout
+ & ensembles,with_dihed_constr,check_conf,histout,constr_dist
double precision Kh(MaxQ,MaxR,MaxT_h,max_parm),
& q0(MaxQ,MaxR,MaxT_h,max_parm),delta,deltrms,deltrgy,fimin,
& f(maxR,maxT_h,max_parm),beta_h(MaxT_h,max_parm)
+ double precision delta_T,startGridT
integer nR(maxT_h,max_parm),snk(MaxR,MaxT_h,max_parm,MaxSlice),
& nT_h(max_parm),maxit,totraj(maxR,max_parm),nRR(maxT_h,max_parm)
+ integer nGridT
logical replica(max_parm),umbrella(max_parm),read_iset(max_parm)
- common /wham/ Kh,q0,f,beta_h,delta,deltrms,deltrgy,fimin,snk,nR,
+ common /wham/ Kh,q0,f,beta_h,delta,deltrms,deltrgy,delta_T,
+ & startGridT,fimin,snk,nR,
& nRR,nT_h,nQ,stot,nparmset,maxit,rescale_mode,replica,umbrella,
- & read_iset,totraj,hamil_rep,separate_parset,iparmprint,myparm
+ & read_iset,totraj,hamil_rep,separate_parset,iparmprint,myparm,
+ & nGridT
parameter (MaxN=100)
integer MaxPrintConf
parameter (MaxPrintConf=1000)
+ integer Max_GridT
+ parameter (Max_GridT=400)
+++ /dev/null
-INSTALL_DIR = /users/software/mpich-1.2.7p1_intel-10.1_em64_ssh
-BIN = ../bin
-FC= ifort
-#OPT = -mcmodel=medium -O3 -ip -w
-OPT = -mcmodel=medium -g -CB
-FFLAGS = ${OPT} -c -I. -I./include_unres -I$(INSTALL_DIR)/include
-LIBS = -L$(INSTALL_DIR)/lib -lmpich -lpmpich xdrf/libxdrf.a
-CPPFLAGS = -DMPI -DLINUX -DUNRES -DSPLITELE -DPROCOR -DPGI -DISNAN -DAMD64 \
- -DCRYST_BOND -DCRYST_THETA -DCRYST_SC
-
-.f.o:
- ${FC} ${FFLAGS} $*.f
-
-.F.o:
- ${FC} ${FFLAGS} ${CPPFLAGS} $*.F
-
-all: make_dbase
-
-objects = \
- wham_multparm.o \
- bxread.o \
- xread.o \
- cxread.o \
- enecalc1.o \
- energy_p_new.o \
- initialize_p.o \
- molread_zs.o \
- openunits.o \
- readrtns.o \
- arcos.o \
- cartder.o \
- cartprint.o \
- chainbuild.o \
- geomout.o \
- icant.o \
- intcor.o \
- int_from_cart.o \
- make_ensemble1.o \
- matmult.o \
- misc.o \
- mygetenv.o \
- parmread.o \
- pinorm.o \
- printmat.o \
- proc_proc.o \
- rescode.o \
- setup_var.o \
- slices.o \
- store_parm.o \
- timing.o \
- wham_calc1.o
-
-objects_compar = \
- readrtns_compar.o \
- readpdb.o fitsq.o contact.o \
- elecont.o contfunc.o cont_frag.o conf_compar.o match_contact.o \
- angnorm.o odlodc.o promienie.o qwolynes.o read_ref_str.o \
- rmscalc.o secondary.o proc_cont.o define_pairs.o mysort.o
-
-make_dbase: ${objects} ${objects_compar}
- cc -o compinfo compinfo.c
- ./compinfo
- ${FC} -c ${FFLAGS} cinfo.f
- $(FC) ${OPT} ${objects} ${objects_compar} cinfo.o \
- ${LIBS} -static-intel -o ${BIN}/wham_multparm-ham_rep-oldparm
-
-clean:
- /bin/rm *.o
--- /dev/null
+Makefile_MPICH_ifort
\ No newline at end of file
--- /dev/null
+INSTALL_DIR = /users/software/mpich-1.2.7p1_intel-10.1_em64_ssh
+BIN = ../../../bin/wham
+FC= ifort
+OPT = -mcmodel=medium -O3 -ip -w
+#OPT = -mcmodel=medium -g -CB
+FFLAGS = ${OPT} -c -I. -I./include_unres -I$(INSTALL_DIR)/include
+LIBS = -L$(INSTALL_DIR)/lib -lmpich -lpmpich xdrf/libxdrf.a
+
+.f.o:
+ ${FC} ${FFLAGS} $*.f
+
+.F.o:
+ ${FC} ${FFLAGS} ${CPPFLAGS} $*.F
+
+objects = \
+ wham_multparm.o \
+ bxread.o \
+ xread.o \
+ cxread.o \
+ enecalc1.o \
+ energy_p_new.o \
+ gnmr1.o \
+ initialize_p.o \
+ molread_zs.o \
+ openunits.o \
+ readrtns.o \
+ arcos.o \
+ cartder.o \
+ cartprint.o \
+ chainbuild.o \
+ geomout.o \
+ icant.o \
+ intcor.o \
+ int_from_cart.o \
+ make_ensemble1.o \
+ matmult.o \
+ misc.o \
+ mygetenv.o \
+ parmread.o \
+ pinorm.o \
+ printmat.o \
+ proc_proc.o \
+ rescode.o \
+ setup_var.o \
+ slices.o \
+ store_parm.o \
+ timing.o \
+ wham_calc1.o
+
+objects_compar = \
+ readrtns_compar.o \
+ readpdb.o fitsq.o contact.o \
+ elecont.o contfunc.o cont_frag.o conf_compar.o match_contact.o \
+ angnorm.o odlodc.o promienie.o qwolynes.o read_ref_str.o \
+ rmscalc.o secondary.o proc_cont.o define_pairs.o mysort.o
+
+GAB: CPPFLAGS = -DMPI -DLINUX -DUNRES -DSPLITELE -DPROCOR -DPGI -DISNAN -DAMD64 \
+ -DCRYST_BOND -DCRYST_THETA -DCRYST_SC
+GAB: ${objects} ${objects_compar} xdrf/libxdrf.a
+ cc -o compinfo compinfo.c
+ ./compinfo
+ ${FC} -c ${FFLAGS} cinfo.f
+ $(FC) ${OPT} ${objects} ${objects_compar} cinfo.o \
+ ${LIBS} -static-intel -o ${BIN}/wham_ifort_MPICH_GAB.exe
+
+E0LL2Y: CPPFLAGS = -DMPI -DLINUX -DUNRES -DSPLITELE -DPROCOR -DPGI -DISNAN -DAMD64
+E0LL2Y: ${objects} ${objects_compar} xdrf/libxdrf.a
+ cc -o compinfo compinfo.c
+ ./compinfo
+ ${FC} -c ${FFLAGS} cinfo.f
+ $(FC) ${OPT} ${objects} ${objects_compar} cinfo.o \
+ ${LIBS} -static-intel -o ${BIN}/wham_ifort_MPICH_E0LL2Y.exe
+
+xdrf/libxdrf.a:
+ cd xdrf && make
+
+
+clean:
+ /bin/rm -f *.o && /bin/rm -f compinfo && cd xdrf && make clean
+
C DO NOT EDIT THIS FILE - IT HAS BEEN GENERATED BY COMPINFO.C
-C 0 0 496
+C 0 0 546
subroutine cinfo
include 'COMMON.IOUNITS'
write(iout,*)'++++ Compile info ++++'
- write(iout,*)'Version 0.0 build 496'
- write(iout,*)'compiled Sun Feb 19 04:44:59 2012'
+ write(iout,*)'Version 0.0 build 546'
+ write(iout,*)'compiled Sun May 20 07:28:19 2012'
write(iout,*)'compiled by adam@matrix.chem.cornell.edu'
write(iout,*)'OS name: Linux '
- write(iout,*)'OS release: 2.6.34.8-68.fc13.x86_64 '
- write(iout,*)'OS version: #1 SMP Thu Feb 17 15:03:58 UTC 2011 '
+ write(iout,*)'OS release: 2.6.34.9-69.fc13.x86_64 '
+ write(iout,*)'OS version: #1 SMP Tue May 3 09:23:03 UTC 2011 '
write(iout,*)'flags:'
write(iout,*)'INSTALL_DIR = /users/software/mpich-1.2.7p1_int...'
- write(iout,*)'BIN = ../bin'
+ write(iout,*)'BIN = ../../../bin/wham'
write(iout,*)'FC= ifort'
- write(iout,*)'OPT = -mcmodel=medium -g -CB'
+ write(iout,*)'OPT = -mcmodel=medium -O3 -ip -w'
write(iout,*)'FFLAGS = ${OPT} -c -I. -I./include_unres -I$(IN...'
write(iout,*)'LIBS = -L$(INSTALL_DIR)/lib -lmpich -lpmpich xd...'
- write(iout,*)'CPPFLAGS = -DMPI -DLINUX -DUNRES -DSPLITELE -DP...'
write(iout,*)'objects = \\'
write(iout,*)' wham_multparm.o \\'
write(iout,*)' bxread.o \\'
write(iout,*)' cxread.o \\'
write(iout,*)' enecalc1.o \\'
write(iout,*)' energy_p_new.o \\'
+ write(iout,*)' gnmr1.o \\'
write(iout,*)' initialize_p.o \\'
write(iout,*)' molread_zs.o \\'
write(iout,*)' openunits.o \\'
write(iout,*)' readrtns_compar.o \\'
write(iout,*)' readpdb.o fitsq.o contact.o \\'
write(iout,*)' elecont.o contfunc.o cont_frag.o conf_c...'
+ write(iout,*)'GAB: CPPFLAGS = -DMPI -DLINUX -DUNRES -DSPLITEL...'
+ write(iout,*)'E0LL2Y: CPPFLAGS = -DMPI -DLINUX -DUNRES -DSPLI...'
write(iout,*)'++++ End of compile info ++++'
return
end
C
implicit real*8 (a-h,o-z)
include 'DIMENSIONS'
- include 'DIMENSIONS.ZSCOPT'
include 'COMMON.SBRIDGE'
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.VAR'
include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
dimension ggg(3)
ehpb=0.0D0
-cd print *,'edis: nhpb=',nhpb,' fbr=',fbr
-cd print *,'link_start=',link_start,' link_end=',link_end
+cd write(iout,*)'edis: nhpb=',nhpb,' fbr=',fbr
+cd write(iout,*)'link_start=',link_start,' link_end=',link_end
if (link_end.eq.0) return
do i=link_start,link_end
C If ihpb(i) and jhpb(i) > NRES, this is a SC-SC distance, otherwise a
iii=ii
jjj=jj
endif
+c write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj,
+c & dhpb(i),dhpb1(i),forcon(i)
C 24/11/03 AL: SS bridges handled separately because of introducing a specific
C distance and angle dependent SS bond potential.
if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then
call ssbond_ene(iii,jjj,eij)
ehpb=ehpb+2*eij
+cd write (iout,*) "eij",eij
+ else if (ii.gt.nres .and. jj.gt.nres) then
+c Restraints from contact prediction
+ dd=dist(ii,jj)
+ if (dhpb1(i).gt.0.0d0) then
+ ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+ fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd
+c write (iout,*) "beta nmr",
+c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+ else
+ dd=dist(ii,jj)
+ rdis=dd-dhpb(i)
+C Get the force constant corresponding to this distance.
+ waga=forcon(i)
+C Calculate the contribution to energy.
+ ehpb=ehpb+waga*rdis*rdis
+c write (iout,*) "beta reg",dd,waga*rdis*rdis
+C
+C Evaluate gradient.
+C
+ fac=waga*rdis/dd
+ endif
+ do j=1,3
+ ggg(j)=fac*(c(j,jj)-c(j,ii))
+ enddo
+ do j=1,3
+ ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
+ ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
+ enddo
+ do k=1,3
+ ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
+ ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
+ enddo
else
C Calculate the distance between the two points and its difference from the
C target distance.
- dd=dist(ii,jj)
- rdis=dd-dhpb(i)
+ dd=dist(ii,jj)
+ if (dhpb1(i).gt.0.0d0) then
+ ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+ fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd
+c write (iout,*) "alph nmr",
+c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+ else
+ rdis=dd-dhpb(i)
C Get the force constant corresponding to this distance.
- waga=forcon(i)
+ waga=forcon(i)
C Calculate the contribution to energy.
- ehpb=ehpb+waga*rdis*rdis
+ ehpb=ehpb+waga*rdis*rdis
+c write (iout,*) "alpha reg",dd,waga*rdis*rdis
C
C Evaluate gradient.
C
- fac=waga*rdis/dd
+ fac=waga*rdis/dd
+ endif
cd print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd,
cd & ' waga=',waga,' fac=',fac
- do j=1,3
- ggg(j)=fac*(c(j,jj)-c(j,ii))
- enddo
+ do j=1,3
+ ggg(j)=fac*(c(j,jj)-c(j,ii))
+ enddo
cd print '(i3,3(1pe14.5))',i,(ggg(j),j=1,3)
C If this is a SC-SC distance, we need to calculate the contributions to the
C Cartesian gradient in the SC vectors (ghpbx).
- if (iii.lt.ii) then
+ if (iii.lt.ii) then
do j=1,3
ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
enddo
- endif
- do j=iii,jjj-1
+ endif
do k=1,3
- ghpbc(k,j)=ghpbc(k,j)+ggg(k)
+ ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
+ ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
enddo
- enddo
endif
enddo
ehpb=0.5D0*ehpb
--- /dev/null
+ double precision function gnmr1(y,ymin,ymax)
+ implicit none
+ double precision y,ymin,ymax
+ double precision wykl /4.0d0/
+ if (y.lt.ymin) then
+ gnmr1=(ymin-y)**wykl/wykl
+ else if (y.gt.ymax) then
+ gnmr1=(y-ymax)**wykl/wykl
+ else
+ gnmr1=0.0d0
+ endif
+ return
+ end
+c------------------------------------------------------------------------------
+ double precision function gnmr1prim(y,ymin,ymax)
+ implicit none
+ double precision y,ymin,ymax
+ double precision wykl /4.0d0/
+ if (y.lt.ymin) then
+ gnmr1prim=-(ymin-y)**(wykl-1)
+ else if (y.gt.ymax) then
+ gnmr1prim=(y-ymax)**(wykl-1)
+ else
+ gnmr1prim=0.0d0
+ endif
+ return
+ end
+c------------------------------------------------------------------------------
+ double precision function harmonic(y,ymax)
+ implicit none
+ double precision y,ymax
+ double precision wykl /2.0d0/
+ harmonic=(y-ymax)**wykl
+ return
+ end
+c-------------------------------------------------------------------------------
+ double precision function harmonicprim(y,ymax)
+ double precision y,ymin,ymax
+ double precision wykl /2.0d0/
+ harmonicprim=(y-ymax)*wykl
+ return
+ end
+c---------------------------------------------------------------------------------
double precision ebr,d0cm,akcm,akth,akct,v1ss,v2ss,v3ss,dhpb,
- & forcon,weidis
- integer ns,nss,nfree,iss,ihpb,jhpb,nhpb,link_start,link_end
+ & dhpb1,forcon,weidis
+ integer ns,nss,nfree,iss,ihpb,jhpb,nhpb,link_start,link_end,
+ & ibecarb
common /sbridge/ ebr,d0cm,akcm,akth,akct,v1ss,v2ss,v3ss,ns,nss,
& nfree,iss(maxss)
- common /links/ dhpb(maxdim),forcon(maxdim),ihpb(maxdim),
- & jhpb(maxdim),nhpb
+ common /links/ dhpb(maxdim),dhpb1(maxdim),forcon(maxdim),
+ & ihpb(maxdim),jhpb(maxdim),ibecarb(maxdim),nhpb
common /restraints/ weidis
common /links_split/ link_start,link_end
if (itype(1).eq.21) nnt=2
if (itype(nres).eq.21) nct=nct-1
write(iout,*) 'NNT=',NNT,' NCT=',NCT
+c Read distance restraints
+ if (constr_dist.gt.0) then
+ call read_dist_constr
+ call hpb_partition
+ endif
+
call setup_var
call init_int_table
if (ns.gt.0) then
return
10 return1
end
+c-------------------------------------------------------------------------------
+ subroutine read_dist_constr
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.CONTROL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.SBRIDGE'
+ integer ifrag_(2,100),ipair_(2,100)
+ double precision wfrag_(100),wpair_(100)
+ character*500 controlcard
+c write (iout,*) "Calling read_dist_constr"
+c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
+c call flush(iout)
+ call card_concat(controlcard,.true.)
+ call readi(controlcard,"NFRAG",nfrag_,0)
+ call readi(controlcard,"NPAIR",npair_,0)
+ call readi(controlcard,"NDIST",ndist_,0)
+ call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
+ call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
+ call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
+ call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
+ call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
+ write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
+ write (iout,*) "IFRAG"
+ do i=1,nfrag_
+ write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
+ enddo
+ write (iout,*) "IPAIR"
+ do i=1,npair_
+ write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
+ enddo
+ call flush(iout)
+ if (.not.refstr .and. nfrag_.gt.0) then
+ write (iout,*)
+ & "ERROR: no reference structure to compute distance restraints"
+ write (iout,*)
+ & "Restraints must be specified explicitly (NDIST=number)"
+ stop
+ endif
+ if (nfrag_.lt.2 .and. npair_.gt.0) then
+ write (iout,*) "ERROR: Less than 2 fragments specified",
+ & " but distance restraints between pairs requested"
+ stop
+ endif
+ call flush(iout)
+ do i=1,nfrag_
+ if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
+ if (ifrag_(2,i).gt.nstart_sup+nsup-1)
+ & ifrag_(2,i)=nstart_sup+nsup-1
+c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
+ call flush(iout)
+ if (wfrag_(i).gt.0.0d0) then
+ do j=ifrag_(1,i),ifrag_(2,i)-1
+ do k=j+1,ifrag_(2,i)
+ write (iout,*) "j",j," k",k
+ ddjk=dist(j,k)
+ if (constr_dist.eq.1) then
+ nhpb=nhpb+1
+ ihpb(nhpb)=j
+ jhpb(nhpb)=k
+ dhpb(nhpb)=ddjk
+ forcon(nhpb)=wfrag_(i)
+ else if (constr_dist.eq.2) then
+ if (ddjk.le.dist_cut) then
+ nhpb=nhpb+1
+ ihpb(nhpb)=j
+ jhpb(nhpb)=k
+ dhpb(nhpb)=ddjk
+ forcon(nhpb)=wfrag_(i)
+ endif
+ else
+ nhpb=nhpb+1
+ ihpb(nhpb)=j
+ jhpb(nhpb)=k
+ dhpb(nhpb)=ddjk
+ forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
+ endif
+ write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
+ & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+ enddo
+ enddo
+ endif
+ enddo
+ do i=1,npair_
+ if (wpair_(i).gt.0.0d0) then
+ ii = ipair_(1,i)
+ jj = ipair_(2,i)
+ if (ii.gt.jj) then
+ itemp=ii
+ ii=jj
+ jj=itemp
+ endif
+ do j=ifrag_(1,ii),ifrag_(2,ii)
+ do k=ifrag_(1,jj),ifrag_(2,jj)
+ nhpb=nhpb+1
+ ihpb(nhpb)=j
+ jhpb(nhpb)=k
+ forcon(nhpb)=wpair_(i)
+ dhpb(nhpb)=dist(j,k)
+ write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
+ & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+ enddo
+ enddo
+ endif
+ enddo
+ do i=1,ndist_
+ read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
+ & ibecarb(i),forcon(nhpb+1)
+ if (forcon(nhpb+1).gt.0.0d0) then
+ nhpb=nhpb+1
+ if (ibecarb(i).gt.0) then
+ ihpb(i)=ihpb(i)+nres
+ jhpb(i)=jhpb(i)+nres
+ endif
+ if (dhpb(nhpb).eq.0.0d0)
+ & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
+ endif
+ enddo
+ do i=1,nhpb
+ write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
+ & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
+ enddo
+ call flush(iout)
+ return
+ end
call reada(controlcard,"DELTA",delta,1.0d-2)
call readi(controlcard,"EINICHECK",einicheck,2)
call reada(controlcard,"DELTRMS",deltrms,5.0d-2)
+ call readi(controlcard,"NGRIDT",NGridT,400)
+ call reada(controlcard,"STARTGRIDT",StartGridT,200.0d0)
+ call reada(controlcard,"DELTA_T",Delta_T,1.0d0)
call reada(controlcard,"DELTRGY",deltrgy,5.0d-2)
call readi(controlcard,"RESCALE",rescale_mode,1)
check_conf=index(controlcard,"NO_CHECK_CONF").eq.0
entfile=index(controlcard,"ENTFILE").gt.0
zscfile=index(controlcard,"ZSCFILE").gt.0
with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
- write (iout,*) "with_dihed_constr ",with_dihed_constr
+ call readi(controlcard,'CONSTR_DIST',constr_dist,0)
+ write (iout,*) "with_dihed_constr ",with_dihed_constr,
+ & " CONSTR_DIST",constr_dist
+ call flush(iout)
return
end
c------------------------------------------------------------------------------
include "DIMENSIONS"
include "DIMENSIONS.ZSCOPT"
include "DIMENSIONS.FREE"
- integer nGridT
- parameter (NGridT=400)
integer MaxBinRms,MaxBinRgy
- parameter (MaxBinRms=1000,MaxBinRgy=1000)
+ parameter (MaxBinRms=100,MaxBinRgy=100)
integer MaxHdim
c parameter (MaxHdim=200000)
parameter (MaxHdim=200)
#ifdef MPI
integer tmax_t,upindE_p
double precision fi_p(MaxR,MaxT_h,Max_Parm)
- double precision sumW_p(0:nGridT,Max_Parm),
- & sumE_p(0:nGridT,Max_Parm),sumEsq_p(0:nGridT,Max_Parm),
- & sumQ_p(MaxQ1,0:nGridT,Max_Parm),
- & sumQsq_p(MaxQ1,0:nGridT,Max_Parm),
- & sumEQ_p(MaxQ1,0:nGridT,Max_Parm),
- & sumEprim_p(MaxQ1,0:nGridT,Max_Parm),
- & sumEbis_p(0:nGridT,Max_Parm)
+ double precision sumW_p(0:Max_GridT,Max_Parm),
+ & sumE_p(0:Max_GridT,Max_Parm),sumEsq_p(0:Max_GridT,Max_Parm),
+ & sumQ_p(MaxQ1,0:Max_GridT,Max_Parm),
+ & sumQsq_p(MaxQ1,0:Max_GridT,Max_Parm),
+ & sumEQ_p(MaxQ1,0:Max_GridT,Max_Parm),
+ & sumEprim_p(MaxQ1,0:Max_GridT,Max_Parm),
+ & sumEbis_p(0:Max_GridT,Max_Parm)
double precision hfin_p(0:MaxHdim,maxT_h),
& hfin_ent_p(0:MaxHdim),histE_p(0:maxindE),sumH,
& hrmsrgy_p(0:MaxBinRgy,0:MaxBinRms,maxT_h)
double precision rgymin_t,rmsmin_t,rgymax_t,rmsmax_t
- double precision potEmin_t,entmin_p,entmax_p
+ double precision potEmin_t_all(maxT_h,Max_Parm),entmin_p,entmax_p
integer histent_p(0:2000)
logical lprint /.true./
#endif
- double precision delta_T /1.0d0/
double precision rgymin,rmsmin,rgymax,rmsmax
double precision sumW(0:NGridT,Max_Parm),sumE(0:NGridT,Max_Parm),
& sumEsq(0:NGridT,Max_Parm),sumQ(MaxQ1,0:NGridT,Max_Parm),
& dd,dd1,dd2,hh,dmin,denom,finorm,avefi,pom,
& hfin(0:MaxHdim,maxT_h),histE(0:maxindE),
& hrmsrgy(0:MaxBinRgy,0:MaxBinRms,maxT_h),
- & potEmin,ent,
+ & potEmin_all(maxT_h,Max_Parm),potEmin,potEmin_min,ent,
& hfin_ent(0:MaxHdim),vmax,aux
double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,
- & eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/,startGridT/200.0d0/,
+ & eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/,
& eplus,eminus,logfac,tanhT,tt
double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
& escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
call flush(iout)
dmin=0.0d0
tmax=0
- potEmin=1.0d10
+ do i=1,nParmset
+ do j=1,nT_h(i)
+ potEmin_all(j,i)=1.0d10
+ enddo
+ enddo
rgymin=1.0d10
rmsmin=1.0d10
rgymax=0.0d0
#else
do i=1,ntot(islice)
#endif
- do j=1,nParmSet
- if (potE(i,j).le.potEmin) potEmin=potE(i,j)
- enddo
if (q(nQ+1,i).lt.rmsmin) rmsmin=q(nQ+1,i)
if (q(nQ+1,i).gt.rmsmax) rmsmax=q(nQ+1,i)
if (q(nQ+2,i).lt.rgymin) rgymin=q(nQ+2,i)
else
ind_point(i)=ind_point(i)+nbin1**(j-1)*ind
endif
-c write (iout,*) "i",i," j",j," q",q(j,i)," ind_point",
-c & ind_point(i)
- call flush(iout)
if (ind_point(i).lt.0 .or. ind_point(i).gt.MaxHdim) then
write (iout,*) "Error - index exceeds range for point",i,
& " q=",q(j,i)," ind",ind_point(i)
call MPI_AllReduce(tmax,tmax_t,1,MPI_INTEGER,MPI_MAX,
& WHAM_COMM,IERROR)
tmax=tmax_t
- call MPI_AllReduce(potEmin,potEmin_t,1,MPI_DOUBLE_PRECISION,
- & MPI_MIN,WHAM_COMM,IERROR)
call MPI_AllReduce(rmsmin,rmsmin_t,1,MPI_DOUBLE_PRECISION,
& MPI_MIN,WHAM_COMM,IERROR)
call MPI_AllReduce(rmsmax,rmsmax_t,1,MPI_DOUBLE_PRECISION,
& MPI_MIN,WHAM_COMM,IERROR)
call MPI_AllReduce(rgymax,rgymax_t,1,MPI_DOUBLE_PRECISION,
& MPI_MAX,WHAM_COMM,IERROR)
- potEmin=potEmin_t/2
rgymin=rgymin_t
rgymax=rgymax_t
rmsmin=rmsmin_t
rmsmax=rmsmax_t
- write (iout,*) "potEmin",potEmin
#endif
rmsmin=deltrms*dint(rmsmin/deltrms)
rmsmax=deltrms*dint(rmsmax/deltrms)
#endif
#ifdef DEBUG
write (iout,*) i,iparm,1.0d0/(beta_h(ib,iparm)*1.987D-3),
- & etot,potEmin
+ & etot
#endif
#ifdef DEBUG
if (iparm.eq.1 .and. ib.eq.1) then
& *(dd-q0(j,kk,ib,iparm))**2
enddo
v(i,kk,ib,iparm)=
- & -beta_h(ib,iparm)*(etot-potEmin+Econstr)
+ & -beta_h(ib,iparm)*(etot+Econstr)
#ifdef DEBUG
write (iout,'(4i5,4e15.5)') i,kk,ib,iparm,
- & etot,potEmin,etot-potEmin,v(i,kk,ib,iparm)
+ & etot,v(i,kk,ib,iparm)
#endif
enddo ! kk
enddo ! ib
#ifdef DEBUG
write (iout,*) "fi before MPI_Reduce me",me,' master',master
do iparm=1,nParmSet
- do ib=1,nT_h(nparmset)
+ do ib=1,nT_h(iparm)
write (iout,*) "iparm",iparm," ib",ib
write (iout,*) "beta=",beta_h(ib,iparm)
write (iout,'(8e15.5)') (fi_p(i,ib,iparm),i=1,nR(ib,iparm))
20 continue
! Now, put together the histograms from all simulations, in order to get the
! unbiased total histogram.
+
+C Determine the minimum free energies
+#ifdef MPI
+ do i=1,scount(me1)
+#else
+ do i=1,ntot(islice)
+#endif
+c write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
+ do iparm=1,nParmSet
+#ifdef DEBUG
+ write (iout,'(2i5,21f8.2)') i,iparm,
+ & (enetb(k,i,iparm),k=1,21)
+#endif
+ call restore_parm(iparm)
+#ifdef DEBUG
+ write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
+ & wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
+ & wtor_d,wsccor,wbond
+#endif
+ do ib=1,nT_h(iparm)
+ if (rescale_mode.eq.1) then
+ quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
+ quotl=1.0d0
+ kfacl=1.0d0
+ do l=1,5
+ quotl1=quotl
+ quotl=quotl*quot
+ kfacl=kfacl*kfac
+ fT(l)=kfacl/(kfacl-1.0d0+quotl)
+ enddo
+#if defined(FUNCTH)
+ tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
+ ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
+#elif defined(FUNCT)
+ ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
+#else
+ ft(6)=1.0d0
+#endif
+ else if (rescale_mode.eq.2) then
+ quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
+ quotl=1.0d0
+ do l=1,5
+ quotl=quotl*quot
+ fT(l)=1.12692801104297249644d0/
+ & dlog(dexp(quotl)+dexp(-quotl))
+ enddo
+#if defined(FUNCTH)
+ tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
+ ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
+#elif defined(FUNCT)
+ ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
+#else
+ ft(6)=1.0d0
+#endif
+c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
+ else if (rescale_mode.eq.0) then
+ do l=1,6
+ fT(l)=1.0d0
+ enddo
+ else
+ write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
+ & rescale_mode
+ call flush(iout)
+ return1
+ endif
+ evdw=enetb(1,i,iparm)
+ evdw_t=enetb(21,i,iparm)
+#ifdef SCP14
+ evdw2_14=enetb(17,i,iparm)
+ evdw2=enetb(2,i,iparm)+evdw2_14
+#else
+ evdw2=enetb(2,i,iparm)
+ evdw2_14=0.0d0
+#endif
+#ifdef SPLITELE
+ ees=enetb(3,i,iparm)
+ evdw1=enetb(16,i,iparm)
+#else
+ ees=enetb(3,i,iparm)
+ evdw1=0.0d0
+#endif
+ ecorr=enetb(4,i,iparm)
+ ecorr5=enetb(5,i,iparm)
+ ecorr6=enetb(6,i,iparm)
+ eel_loc=enetb(7,i,iparm)
+ eello_turn3=enetb(8,i,iparm)
+ eello_turn4=enetb(9,i,iparm)
+ eturn6=enetb(10,i,iparm)
+ ebe=enetb(11,i,iparm)
+ escloc=enetb(12,i,iparm)
+ etors=enetb(13,i,iparm)
+ etors_d=enetb(14,i,iparm)
+ ehpb=enetb(15,i,iparm)
+ estr=enetb(18,i,iparm)
+ esccor=enetb(19,i,iparm)
+ edihcnstr=enetb(20,i,iparm)
+#ifdef DEBUG
+ write (iout,'(3i5,6f5.2,14f12.3)') i,ib,iparm,(ft(l),l=1,6),
+ & evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,
+ & etors,etors_d,eello_turn3,eello_turn4,esccor,edihcnstr
+#endif
+
+#ifdef SPLITELE
+ etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
+ & +wvdwpp*evdw1
+ & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
+ & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
+ & +ft(2)*wturn3*eello_turn3
+ & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
+ & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
+ & +wbond*estr
+#else
+ etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
+ & +ft(1)*welec*(ees+evdw1)
+ & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
+ & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
+ & +ft(2)*wturn3*eello_turn3
+ & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
+ & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
+ & +wbond*estr
+#endif
+c write (iout,*) "i",i," ib",ib,
+c & " temp",1.0d0/(1.987d-3*beta_h(ib,iparm))," etot",etot,
+c & " entfac",entfac(i)
+ etot=etot-entfac(i)/beta_h(ib,iparm)
+ if(etot.lt.potEmin_all(ib,iparm)) potEmin_all(ib,iparm)=etot
+c write (iout,*) "efree",etot," potEmin",potEmin_all(ib,iparm)
+ enddo ! ib
+ enddo ! iparm
+ enddo ! i
+#ifdef DEBUG
+ write (iout,*) "The potEmin array before reduction"
+ do i=1,nParmSet
+ write (iout,*) "Parameter set",i
+ do j=1,nT_h(i)
+ write (iout,*) j,PotEmin_all(j,i)
+ enddo
+ enddo
+ write (iout,*) "potEmin_min",potEmin_min
+#endif
+#ifdef MPI
+C Determine the minimum energes for all parameter sets and temperatures
+ call MPI_AllReduce(potEmin_all(1,1),potEmin_t_all(1,1),
+ & maxT_h*nParmSet,MPI_DOUBLE_PRECISION,MPI_MIN,WHAM_COMM,IERROR)
+ do i=1,nParmSet
+ do j=1,nT_h(i)
+ potEmin_all(j,i)=potEmin_t_all(j,i)
+ enddo
+ enddo
+#endif
+ potEmin_min=potEmin_all(1,1)
+ do i=1,nParmSet
+ do j=1,nT_h(i)
+ if (potEmin_all(j,i).lt.potEmin_min)
+ & potEmin_min=potEmin_all(j,i)
+ enddo
+ enddo
+#ifdef DEBUG
+ write (iout,*) "The potEmin array"
+ do i=1,nParmSet
+ write (iout,*) "Parameter set",i
+ do j=1,nT_h(i)
+ write (iout,*) j,PotEmin_all(j,i)
+ enddo
+ enddo
+ write (iout,*) "potEmin_min",potEmin_min
+#endif
+#undef DEBUG
+
#ifdef MPI
do t=0,tmax
hfin_ent_p(t)=0.0d0
#else
hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
#endif
-c write (iout,'(2i5,20f8.2)') t,t,(enetb(k,t,iparm),k=1,18)
call restore_parm(iparm)
evdw=enetb(21,t,iparm)
evdw_t=enetb(1,t,iparm)
estr=enetb(18,t,iparm)
esccor=enetb(19,t,iparm)
edihcnstr=enetb(20,t,iparm)
- edihcnstr=0.0d0
do k=0,nGridT
betaT=startGridT+k*delta_T
temper=betaT
c write (iout,*) "ftprim",ftprim
c write (iout,*) "ftbis",ftbis
betaT=1.0d0/(1.987D-3*betaT)
+ if (betaT.ge.beta_h(1,iparm)) then
+ potEmin=potEmin_all(1,iparm)
+c write(iout,*) "first",temper,potEmin
+ else if (betaT.lt.beta_h(nT_h(iparm),iparm)) then
+ potEmin=potEmin_all(nT_h(iparm),iparm)
+c write (iout,*) "last",temper,potEmin
+ else
+ do l=1,nT_h(iparm)-1
+ if (betaT.le.beta_h(l,iparm) .and.
+ & betaT.gt.beta_h(l+1,iparm)) then
+ potEmin=potEmin_all(l,iparm)
+c write (iout,*) "l",l,
+c & betaT,1.0d0/(1.987D-3*beta_h(l,iparm)),
+c & 1.0d0/(1.987D-3*beta_h(l+1,iparm)),temper,potEmin
+ exit
+ endif
+ enddo
+ endif
+c write (iout,*) ib," PotEmin",potEmin
#ifdef SPLITELE
etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
& +wvdwpp*evdw1
#endif
weight=dexp(-betaT*(etot-potEmin)+entfac(t))
#ifdef DEBUG
- write (iout,*) "iparm",iparm," t",t," betaT",betaT,
+ write (iout,*) "iparm",iparm," t",t," temper",temper,
& " etot",etot," entfac",entfac(t),
+ & " efree",etot-entfac(t)/betaT," potEmin",potEmin,
+ & " boltz",-betaT*(etot-potEmin)+entfac(t),
& " weight",weight," ebis",ebis
#endif
etot=etot-temper*eprim
endif
#ifdef MPI
do ib=1,nT_h(iparm)
+ potEmin=potEmin_all(ib,iparm)
expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
hfin_p(ind,ib)=hfin_p(ind,ib)+
& dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
enddo
#else
do ib=1,nT_h(iparm)
+ potEmin=potEmin_all(ib,iparm)
expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
hfin(ind,ib)=hfin(ind,ib)+
& dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
write (iout,'(a,i3)') "Parameter set",iparm
endif
do i=0,NGridT
+ betaT=1.0d0/(1.987D-3*(startGridT+i*delta_T))
+ if (betaT.ge.beta_h(1,iparm)) then
+ potEmin=potEmin_all(1,iparm)
+ else if (betaT.lt.beta_h(nT_h(iparm),iparm)) then
+ potEmin=potEmin_all(nT_h(iparm),iparm)
+ else
+ do l=1,nT_h(iparm)-1
+ if (betaT.le.beta_h(l,iparm) .and.
+ & betaT.gt.beta_h(l+1,iparm)) then
+ potEmin=potEmin_all(l,iparm)
+ exit
+ endif
+ enddo
+ endif
sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/
& sumW(i,iparm)
--- /dev/null
+ subroutine WHAM_CALC(islice,*)
+! Weighed Histogram Analysis Method (WHAM) code
+! Written by A. Liwo based on the work of Kumar et al.,
+! J.Comput.Chem., 13, 1011 (1992)
+!
+! 2/1/05 Multiple temperatures allowed.
+! 2/2/05 Free energies calculated directly from data points
+! acc. to Eq. (21) of Kumar et al.; final histograms also
+! constructed based on this equation.
+! 2/12/05 Multiple parameter sets included
+!
+! 2/2/05 Parallel version
+ implicit none
+ include "DIMENSIONS"
+ include "DIMENSIONS.ZSCOPT"
+ include "DIMENSIONS.FREE"
+ integer nGridT
+ parameter (NGridT=400)
+ integer MaxBinRms,MaxBinRgy
+ parameter (MaxBinRms=100,MaxBinRgy=100)
+ integer MaxHdim
+c parameter (MaxHdim=200000)
+ parameter (MaxHdim=200)
+ integer maxinde
+ parameter (maxinde=200)
+#ifdef MPI
+ include "mpif.h"
+ include "COMMON.MPI"
+ integer ierror,errcode,status(MPI_STATUS_SIZE)
+#endif
+ include "COMMON.CONTROL"
+ include "COMMON.IOUNITS"
+ include "COMMON.FREE"
+ include "COMMON.ENERGIES"
+ include "COMMON.FFIELD"
+ include "COMMON.SBRIDGE"
+ include "COMMON.PROT"
+ include "COMMON.ENEPS"
+ integer MaxPoint,MaxPointProc
+ parameter (MaxPoint=MaxStr,
+ & MaxPointProc=MaxStr_Proc)
+ double precision finorm_max,potfac,entmin,entmax,expfac,vf
+ parameter (finorm_max=1.0d0)
+ integer islice
+ integer i,ii,j,jj,k,kk,l,m,ind,iter,t,tmax,ient,ientmax,iln
+ integer start,end,iharm,ib,iib,nbin1,nbin,nbin_rms,nbin_rgy,
+ & nbin_rmsrgy,liczba,iparm,nFi,indrgy,indrms
+ integer htot(0:MaxHdim),histent(0:2000)
+ double precision v(MaxPointProc,MaxR,MaxT_h,Max_Parm)
+ double precision energia(0:max_ene)
+#ifdef MPI
+ integer tmax_t,upindE_p
+ double precision fi_p(MaxR,MaxT_h,Max_Parm)
+ double precision sumW_p(0:nGridT,Max_Parm),
+ & sumE_p(0:nGridT,Max_Parm),sumEsq_p(0:nGridT,Max_Parm),
+ & sumQ_p(MaxQ1,0:nGridT,Max_Parm),
+ & sumQsq_p(MaxQ1,0:nGridT,Max_Parm),
+ & sumEQ_p(MaxQ1,0:nGridT,Max_Parm),
+ & sumEprim_p(MaxQ1,0:nGridT,Max_Parm),
+ & sumEbis_p(0:nGridT,Max_Parm)
+ double precision hfin_p(0:MaxHdim,maxT_h),
+ & hfin_ent_p(0:MaxHdim),histE_p(0:maxindE),sumH,
+ & hrmsrgy_p(0:MaxBinRgy,0:MaxBinRms,maxT_h)
+ double precision rgymin_t,rmsmin_t,rgymax_t,rmsmax_t
+ double precision potEmin_t,entmin_p,entmax_p
+ integer histent_p(0:2000)
+ logical lprint /.true./
+#endif
+ double precision delta_T /1.0d0/
+ double precision rgymin,rmsmin,rgymax,rmsmax
+ double precision sumW(0:NGridT,Max_Parm),sumE(0:NGridT,Max_Parm),
+ & sumEsq(0:NGridT,Max_Parm),sumQ(MaxQ1,0:NGridT,Max_Parm),
+ & sumQsq(MaxQ1,0:NGridT,Max_Parm),sumEQ(MaxQ1,0:NGridT,Max_Parm),
+ & sumEprim(0:NGridT,Max_Parm),sumEbis(0:NGridT,Max_Parm),betaT,
+ & weight,econstr
+ double precision fi(MaxR,maxT_h,Max_Parm),
+ & dd,dd1,dd2,hh,dmin,denom,finorm,avefi,pom,
+ & hfin(0:MaxHdim,maxT_h),histE(0:maxindE),
+ & hrmsrgy(0:MaxBinRgy,0:MaxBinRms,maxT_h),
+ & potEmin,ent,
+ & hfin_ent(0:MaxHdim),vmax,aux
+ double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,
+ & eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/,startGridT/200.0d0/,
+ & eplus,eminus,logfac,tanhT,tt
+ double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
+ & escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
+ & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor
+
+ integer ind_point(maxpoint),upindE,indE
+ character*16 plik
+ character*1 licz1
+ character*2 licz2
+ character*3 licz3
+ character*128 nazwa
+ integer ilen
+ external ilen
+
+ write(licz2,'(bz,i2.2)') islice
+ nbin1 = 1.0d0/delta
+ write (iout,'(//80(1h-)/"Solving WHAM equations for slice",
+ & i2/80(1h-)//)') islice
+ write (iout,*) "delta",delta," nbin1",nbin1
+ write (iout,*) "MaxN",MaxN," MaxQ",MaxQ," MaHdim",MaxHdim
+ call flush(iout)
+ dmin=0.0d0
+ tmax=0
+ potEmin=1.0d10
+ rgymin=1.0d10
+ rmsmin=1.0d10
+ rgymax=0.0d0
+ rmsmax=0.0d0
+ do t=0,MaxN
+ htot(t)=0
+ enddo
+#ifdef MPI
+ do i=1,scount(me1)
+#else
+ do i=1,ntot(islice)
+#endif
+ do j=1,nParmSet
+ if (potE(i,j).le.potEmin) potEmin=potE(i,j)
+ enddo
+ if (q(nQ+1,i).lt.rmsmin) rmsmin=q(nQ+1,i)
+ if (q(nQ+1,i).gt.rmsmax) rmsmax=q(nQ+1,i)
+ if (q(nQ+2,i).lt.rgymin) rgymin=q(nQ+2,i)
+ if (q(nQ+2,i).gt.rgymax) rgymax=q(nQ+2,i)
+ ind_point(i)=0
+ do j=nQ,1,-1
+ ind=(q(j,i)-dmin+1.0d-8)/delta
+ if (j.eq.1) then
+ ind_point(i)=ind_point(i)+ind
+ else
+ ind_point(i)=ind_point(i)+nbin1**(j-1)*ind
+ endif
+c write (iout,*) "i",i," j",j," q",q(j,i)," ind_point",
+c & ind_point(i)
+ call flush(iout)
+ if (ind_point(i).lt.0 .or. ind_point(i).gt.MaxHdim) then
+ write (iout,*) "Error - index exceeds range for point",i,
+ & " q=",q(j,i)," ind",ind_point(i)
+#ifdef MPI
+ write (iout,*) "Processor",me1
+ call flush(iout)
+ call MPI_Abort(MPI_COMM_WORLD, Ierror, Errcode )
+#endif
+ stop
+ endif
+ enddo ! j
+ if (ind_point(i).gt.tmax) tmax=ind_point(i)
+ htot(ind_point(i))=htot(ind_point(i))+1
+#ifdef DEBUG
+ write (iout,*) "i",i,"q",(q(j,i),j=1,nQ)," ind",ind_point(i),
+ & " htot",htot(ind_point(i))
+ call flush(iout)
+#endif
+ enddo ! i
+ call flush(iout)
+
+ nbin=nbin1**nQ-1
+ write (iout,'(a)') "Numbers of counts in Q bins"
+ do t=0,tmax
+ if (htot(t).gt.0) then
+ write (iout,'(i15,$)') t
+ liczba=t
+ do j=1,nQ
+ jj = mod(liczba,nbin1)
+ liczba=liczba/nbin1
+ write (iout,'(i5,$)') jj
+ enddo
+ write (iout,'(i8)') htot(t)
+ endif
+ enddo
+ do iparm=1,nParmSet
+ write (iout,'(a,i3)') "Number of data points for parameter set",
+ & iparm
+ write (iout,'(i7,$)') ((snk(m,ib,iparm,islice),m=1,nr(ib,iparm)),
+ & ib=1,nT_h(iparm))
+ write (iout,'(i8)') stot(islice)
+ write (iout,'(a)')
+ enddo
+ call flush(iout)
+
+#ifdef MPI
+ call MPI_AllReduce(tmax,tmax_t,1,MPI_INTEGER,MPI_MAX,
+ & WHAM_COMM,IERROR)
+ tmax=tmax_t
+ call MPI_AllReduce(potEmin,potEmin_t,1,MPI_DOUBLE_PRECISION,
+ & MPI_MIN,WHAM_COMM,IERROR)
+ call MPI_AllReduce(rmsmin,rmsmin_t,1,MPI_DOUBLE_PRECISION,
+ & MPI_MIN,WHAM_COMM,IERROR)
+ call MPI_AllReduce(rmsmax,rmsmax_t,1,MPI_DOUBLE_PRECISION,
+ & MPI_MAX,WHAM_COMM,IERROR)
+ call MPI_AllReduce(rgymin,rgymin_t,1,MPI_DOUBLE_PRECISION,
+ & MPI_MIN,WHAM_COMM,IERROR)
+ call MPI_AllReduce(rgymax,rgymax_t,1,MPI_DOUBLE_PRECISION,
+ & MPI_MAX,WHAM_COMM,IERROR)
+ potEmin=potEmin_t/2
+ rgymin=rgymin_t
+ rgymax=rgymax_t
+ rmsmin=rmsmin_t
+ rmsmax=rmsmax_t
+ write (iout,*) "potEmin",potEmin
+#endif
+ rmsmin=deltrms*dint(rmsmin/deltrms)
+ rmsmax=deltrms*dint(rmsmax/deltrms)
+ rgymin=deltrms*dint(rgymin/deltrgy)
+ rgymax=deltrms*dint(rgymax/deltrgy)
+ nbin_rms=(rmsmax-rmsmin)/deltrms
+ nbin_rgy=(rgymax-rgymin)/deltrgy
+ write (iout,*) "rmsmin",rmsmin," rmsmax",rmsmax," rgymin",rgymin,
+ & " rgymax",rgymax," nbin_rms",nbin_rms," nbin_rgy",nbin_rgy
+ nFi=0
+ do i=1,nParmSet
+ do j=1,nT_h(i)
+ nFi=nFi+nR(j,i)
+ enddo
+ enddo
+ write (iout,*) "nFi",nFi
+! Compute the Boltzmann factor corresponing to restrain potentials in different
+! simulations.
+#ifdef MPI
+ do i=1,scount(me1)
+#else
+ do i=1,ntot(islice)
+#endif
+c write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
+ do iparm=1,nParmSet
+#ifdef DEBUG
+ write (iout,'(2i5,21f8.2)') i,iparm,
+ & (enetb(k,i,iparm),k=1,21)
+#endif
+ call restore_parm(iparm)
+#ifdef DEBUG
+ write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
+ & wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
+ & wtor_d,wsccor,wbond
+#endif
+ do ib=1,nT_h(iparm)
+ if (rescale_mode.eq.1) then
+ quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
+ quotl=1.0d0
+ kfacl=1.0d0
+ do l=1,5
+ quotl1=quotl
+ quotl=quotl*quot
+ kfacl=kfacl*kfac
+ fT(l)=kfacl/(kfacl-1.0d0+quotl)
+ enddo
+#if defined(FUNCTH)
+ tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
+ ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
+#elif defined(FUNCT)
+ ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
+#else
+ ft(6)=1.0d0
+#endif
+ else if (rescale_mode.eq.2) then
+ quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
+ quotl=1.0d0
+ do l=1,5
+ quotl=quotl*quot
+ fT(l)=1.12692801104297249644d0/
+ & dlog(dexp(quotl)+dexp(-quotl))
+ enddo
+#if defined(FUNCTH)
+ tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
+ ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
+#elif defined(FUNCT)
+ ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
+#else
+ ft(6)=1.0d0
+#endif
+c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
+ else if (rescale_mode.eq.0) then
+ do l=1,6
+ fT(l)=1.0d0
+ enddo
+ else
+ write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
+ & rescale_mode
+ call flush(iout)
+ return1
+ endif
+ evdw=enetb(1,i,iparm)
+ evdw_t=enetb(21,i,iparm)
+#ifdef SCP14
+ evdw2_14=enetb(17,i,iparm)
+ evdw2=enetb(2,i,iparm)+evdw2_14
+#else
+ evdw2=enetb(2,i,iparm)
+ evdw2_14=0.0d0
+#endif
+#ifdef SPLITELE
+ ees=enetb(3,i,iparm)
+ evdw1=enetb(16,i,iparm)
+#else
+ ees=enetb(3,i,iparm)
+ evdw1=0.0d0
+#endif
+ ecorr=enetb(4,i,iparm)
+ ecorr5=enetb(5,i,iparm)
+ ecorr6=enetb(6,i,iparm)
+ eel_loc=enetb(7,i,iparm)
+ eello_turn3=enetb(8,i,iparm)
+ eello_turn4=enetb(9,i,iparm)
+ eturn6=enetb(10,i,iparm)
+ ebe=enetb(11,i,iparm)
+ escloc=enetb(12,i,iparm)
+ etors=enetb(13,i,iparm)
+ etors_d=enetb(14,i,iparm)
+ ehpb=enetb(15,i,iparm)
+ estr=enetb(18,i,iparm)
+ esccor=enetb(19,i,iparm)
+ edihcnstr=enetb(20,i,iparm)
+#ifdef DEBUG
+ write (iout,'(3i5,6f5.2,14f12.3)') i,ib,iparm,(ft(l),l=1,6),
+ & evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,
+ & etors,etors_d,eello_turn3,eello_turn4,esccor
+#endif
+
+#ifdef SPLITELE
+ etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
+ & +wvdwpp*evdw1
+ & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
+ & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
+ & +ft(2)*wturn3*eello_turn3
+ & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
+ & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
+ & +wbond*estr
+#else
+ etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
+ & +ft(1)*welec*(ees+evdw1)
+ & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
+ & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
+ & +ft(2)*wturn3*eello_turn3
+ & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
+ & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
+ & +wbond*estr
+#endif
+#ifdef DEBUG
+ write (iout,*) i,iparm,1.0d0/(beta_h(ib,iparm)*1.987D-3),
+ & etot,potEmin
+#endif
+#ifdef DEBUG
+ if (iparm.eq.1 .and. ib.eq.1) then
+ write (iout,*)"Conformation",i
+ energia(0)=etot
+ do k=1,max_ene
+ energia(k)=enetb(k,i,iparm)
+ enddo
+ call enerprint(energia(0),fT)
+ endif
+#endif
+ do kk=1,nR(ib,iparm)
+ Econstr=0.0d0
+ do j=1,nQ
+ dd = q(j,i)
+ Econstr=Econstr+Kh(j,kk,ib,iparm)
+ & *(dd-q0(j,kk,ib,iparm))**2
+ enddo
+ v(i,kk,ib,iparm)=
+ & -beta_h(ib,iparm)*(etot-potEmin+Econstr)
+#ifdef DEBUG
+ write (iout,'(4i5,4e15.5)') i,kk,ib,iparm,
+ & etot,potEmin,etot-potEmin,v(i,kk,ib,iparm)
+#endif
+ enddo ! kk
+ enddo ! ib
+ enddo ! iparm
+ enddo ! i
+! Simple iteration to calculate free energies corresponding to all simulation
+! runs.
+ do iter=1,maxit
+
+! Compute new free-energy values corresponding to the righ-hand side of the
+! equation and their derivatives.
+ write (iout,*) "------------------------fi"
+#ifdef MPI
+ do t=1,scount(me1)
+#else
+ do t=1,ntot(islice)
+#endif
+ vmax=-1.0d+20
+ do i=1,nParmSet
+ do k=1,nT_h(i)
+ do l=1,nR(k,i)
+ vf=v(t,l,k,i)+f(l,k,i)
+ if (vf.gt.vmax) vmax=vf
+ enddo
+ enddo
+ enddo
+ denom=0.0d0
+ do i=1,nParmSet
+ do k=1,nT_h(i)
+ do l=1,nR(k,i)
+ aux=f(l,k,i)+v(t,l,k,i)-vmax
+ if (aux.gt.-200.0d0)
+ & denom=denom+snk(l,k,i,islice)*dexp(aux)
+ enddo
+ enddo
+ enddo
+ entfac(t)=-dlog(denom)-vmax
+#ifdef DEBUG
+ write (iout,*) t,"vmax",vmax," denom",denom,"entfac",entfac(t)
+#endif
+ enddo
+ do iparm=1,nParmSet
+ do iib=1,nT_h(iparm)
+ do ii=1,nR(iib,iparm)
+#ifdef MPI
+ fi_p(ii,iib,iparm)=0.0d0
+ do t=1,scount(me)
+ fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm)
+ & +dexp(v(t,ii,iib,iparm)+entfac(t))
+#ifdef DEBUG
+ write (iout,'(4i5,3e15.5)') t,ii,iib,iparm,
+ & v(t,ii,iib,iparm),entfac(t),fi_p(ii,iib,iparm)
+#endif
+ enddo
+#else
+ fi(ii,iib,iparm)=0.0d0
+ do t=1,ntot(islice)
+ fi(ii,iib,iparm)=fi(ii,iib,iparm)
+ & +dexp(v(t,ii,iib,iparm)+entfac(t))
+ enddo
+#endif
+ enddo ! ii
+ enddo ! iib
+ enddo ! iparm
+
+#ifdef MPI
+#ifdef DEBUG
+ write (iout,*) "fi before MPI_Reduce me",me,' master',master
+ do iparm=1,nParmSet
+ do ib=1,nT_h(nparmset)
+ write (iout,*) "iparm",iparm," ib",ib
+ write (iout,*) "beta=",beta_h(ib,iparm)
+ write (iout,'(8e15.5)') (fi_p(i,ib,iparm),i=1,nR(ib,iparm))
+ enddo
+ enddo
+#endif
+ write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet,
+ & maxR*MaxT_h*nParmSet
+ write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD,
+ & " WHAM_COMM",WHAM_COMM
+ call MPI_Reduce(fi_p(1,1,1),fi(1,1,1),maxR*MaxT_h*nParmSet,
+ & MPI_DOUBLE_PRECISION,
+ & MPI_SUM,Master,WHAM_COMM,IERROR)
+#ifdef DEBUG
+ write (iout,*) "fi after MPI_Reduce nparmset",nparmset
+ do iparm=1,nParmSet
+ write (iout,*) "iparm",iparm
+ do ib=1,nT_h(iparm)
+ write (iout,*) "beta=",beta_h(ib,iparm)
+ write (iout,'(8e15.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
+ enddo
+ enddo
+#endif
+ if (me1.eq.Master) then
+#endif
+ avefi=0.0d0
+ do iparm=1,nParmSet
+ do ib=1,nT_h(iparm)
+ do i=1,nR(ib,iparm)
+ fi(i,ib,iparm)=-dlog(fi(i,ib,iparm))
+ avefi=avefi+fi(i,ib,iparm)
+ enddo
+ enddo
+ enddo
+ avefi=avefi/nFi
+ do iparm=1,nParmSet
+ write (iout,*) "Parameter set",iparm
+ do ib =1,nT_h(iparm)
+ write (iout,*) "beta=",beta_h(ib,iparm)
+ do i=1,nR(ib,iparm)
+ fi(i,ib,iparm)=fi(i,ib,iparm)-avefi
+ enddo
+ write (iout,'(8f10.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
+ write (iout,'(8f10.5)') (f(i,ib,iparm),i=1,nR(ib,iparm))
+ enddo
+ enddo
+
+! Compute the norm of free-energy increments.
+ finorm=0.0d0
+ do iparm=1,nParmSet
+ do ib=1,nT_h(iparm)
+ do i=1,nR(ib,iparm)
+ finorm=finorm+dabs(fi(i,ib,iparm)-f(i,ib,iparm))
+ f(i,ib,iparm)=fi(i,ib,iparm)
+ enddo
+ enddo
+ enddo
+
+ write (iout,*) 'Iteration',iter,' finorm',finorm
+
+#ifdef MPI
+ endif
+ call MPI_Bcast(f(1,1,1),MaxR*MaxT_h*nParmSet,
+ & MPI_DOUBLE_PRECISION,Master,
+ & WHAM_COMM,IERROR)
+ call MPI_Bcast(finorm,1,MPI_DOUBLE_PRECISION,Master,
+ & WHAM_COMM,IERROR)
+#endif
+! Exit, if the increment norm is smaller than pre-assigned tolerance.
+ if (finorm.lt.fimin) then
+ write (iout,*) 'Iteration converged'
+ goto 20
+ endif
+
+ enddo ! iter
+
+ 20 continue
+! Now, put together the histograms from all simulations, in order to get the
+! unbiased total histogram.
+#ifdef MPI
+ do t=0,tmax
+ hfin_ent_p(t)=0.0d0
+ enddo
+#else
+ do t=0,tmax
+ hfin_ent(t)=0.0d0
+ enddo
+#endif
+ write (iout,*) "--------------hist"
+#ifdef MPI
+ do iparm=1,nParmSet
+ do i=0,nGridT
+ sumW_p(i,iparm)=0.0d0
+ sumE_p(i,iparm)=0.0d0
+ sumEbis_p(i,iparm)=0.0d0
+ sumEsq_p(i,iparm)=0.0d0
+ do j=1,nQ+2
+ sumQ_p(j,i,iparm)=0.0d0
+ sumQsq_p(j,i,iparm)=0.0d0
+ sumEQ_p(j,i,iparm)=0.0d0
+ enddo
+ enddo
+ enddo
+ upindE_p=0
+#else
+ do iparm=1,nParmSet
+ do i=0,nGridT
+ sumW(i,iparm)=0.0d0
+ sumE(i,iparm)=0.0d0
+ sumEbis(i,iparm)=0.0d0
+ sumEsq(i,iparm)=0.0d0
+ do j=1,nQ+2
+ sumQ(j,i,iparm)=0.0d0
+ sumQsq(j,i,iparm)=0.0d0
+ sumEQ(j,i,iparm)=0.0d0
+ enddo
+ enddo
+ enddo
+ upindE=0
+#endif
+c 8/26/05 entropy distribution
+#ifdef MPI
+ entmin_p=1.0d10
+ entmax_p=-1.0d10
+ do t=1,scount(me1)
+c ent=-dlog(entfac(t))
+ ent=entfac(t)
+ if (ent.lt.entmin_p) entmin_p=ent
+ if (ent.gt.entmax_p) entmax_p=ent
+ enddo
+ write (iout,*) "entmin",entmin_p," entmax",entmax_p
+ call flush(iout)
+ call MPI_Allreduce(entmin_p,entmin,1,MPI_DOUBLE_PRECISION,MPI_MIN,
+ & WHAM_COMM,IERROR)
+ call MPI_Allreduce(entmax_p,entmax,1,MPI_DOUBLE_PRECISION,MPI_MAX,
+ & WHAM_COMM,IERROR)
+ ientmax=entmax-entmin
+ if (ientmax.gt.2000) ientmax=2000
+ write (iout,*) "entmin",entmin," entmax",entmax," ientmax",ientmax
+ call flush(iout)
+ do t=1,scount(me1)
+c ient=-dlog(entfac(t))-entmin
+ ient=entfac(t)-entmin
+ if (ient.le.2000) histent_p(ient)=histent_p(ient)+1
+ enddo
+ call MPI_Allreduce(histent_p(0),histent(0),ientmax+1,MPI_INTEGER,
+ & MPI_SUM,WHAM_COMM,IERROR)
+ if (me1.eq.Master) then
+ write (iout,*) "Entropy histogram"
+ do i=0,ientmax
+ write(iout,'(f15.4,i10)') entmin+i,histent(i)
+ enddo
+ endif
+#else
+ entmin=1.0d10
+ entmax=-1.0d10
+ do t=1,ntot(islice)
+ ent=entfac(t)
+ if (ent.lt.entmin) entmin=ent
+ if (ent.gt.entmax) entmax=ent
+ enddo
+ ientmax=-dlog(entmax)-entmin
+ if (ientmax.gt.2000) ientmax=2000
+ do t=1,ntot(islice)
+ ient=entfac(t)-entmin
+ if (ient.le.2000) histent(ient)=histent(ient)+1
+ enddo
+ write (iout,*) "Entropy histogram"
+ do i=0,ientmax
+ write(iout,'(2f15.4)') entmin+i,histent(i)
+ enddo
+#endif
+
+#ifdef MPI
+c write (iout,*) "me1",me1," scount",scount(me1)
+
+ do iparm=1,nParmSet
+
+#ifdef MPI
+ do ib=1,nT_h(iparm)
+ do t=0,tmax
+ hfin_p(t,ib)=0.0d0
+ enddo
+ enddo
+ do i=1,maxindE
+ histE_p(i)=0.0d0
+ enddo
+#else
+ do ib=1,nT_h(iparm)
+ do t=0,tmax
+ hfin(t,ib)=0.0d0
+ enddo
+ enddo
+ do i=1,maxindE
+ histE(i)=0.0d0
+ enddo
+#endif
+ do ib=1,nT_h(iparm)
+ do i=0,MaxBinRms
+ do j=0,MaxBinRgy
+ hrmsrgy(j,i,ib)=0.0d0
+#ifdef MPI
+ hrmsrgy_p(j,i,ib)=0.0d0
+#endif
+ enddo
+ enddo
+ enddo
+
+ do t=1,scount(me1)
+#else
+ do t=1,ntot(islice)
+#endif
+ ind = ind_point(t)
+#ifdef MPI
+ hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
+#else
+ hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
+#endif
+c write (iout,'(2i5,20f8.2)') t,t,(enetb(k,t,iparm),k=1,18)
+ call restore_parm(iparm)
+ evdw=enetb(21,t,iparm)
+ evdw_t=enetb(1,t,iparm)
+#ifdef SCP14
+ evdw2_14=enetb(17,t,iparm)
+ evdw2=enetb(2,t,iparm)+evdw2_14
+#else
+ evdw2=enetb(2,t,iparm)
+ evdw2_14=0.0d0
+#endif
+#ifdef SPLITELE
+ ees=enetb(3,t,iparm)
+ evdw1=enetb(16,t,iparm)
+#else
+ ees=enetb(3,t,iparm)
+ evdw1=0.0d0
+#endif
+ ecorr=enetb(4,t,iparm)
+ ecorr5=enetb(5,t,iparm)
+ ecorr6=enetb(6,t,iparm)
+ eel_loc=enetb(7,t,iparm)
+ eello_turn3=enetb(8,t,iparm)
+ eello_turn4=enetb(9,t,iparm)
+ eturn6=enetb(10,t,iparm)
+ ebe=enetb(11,t,iparm)
+ escloc=enetb(12,t,iparm)
+ etors=enetb(13,t,iparm)
+ etors_d=enetb(14,t,iparm)
+ ehpb=enetb(15,t,iparm)
+ estr=enetb(18,t,iparm)
+ esccor=enetb(19,t,iparm)
+ edihcnstr=enetb(20,t,iparm)
+ edihcnstr=0.0d0
+ do k=0,nGridT
+ betaT=startGridT+k*delta_T
+ temper=betaT
+c fT=T0/betaT
+c ft=2*T0/(T0+betaT)
+ if (rescale_mode.eq.1) then
+ quot=betaT/T0
+ quotl=1.0d0
+ kfacl=1.0d0
+ do l=1,5
+ quotl1=quotl
+ quotl=quotl*quot
+ kfacl=kfacl*kfac
+ denom=kfacl-1.0d0+quotl
+ fT(l)=kfacl/denom
+ ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
+ ftbis(l)=l*kfacl*quotl1*
+ & (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
+ enddo
+#if defined(FUNCTH)
+ ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
+ & 320.0d0
+ ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
+ ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
+ & /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
+#elif defined(FUNCT)
+ fT(6)=betaT/T0
+ ftprim(6)=1.0d0/T0
+ ftbis(6)=0.0d0
+#else
+ fT(6)=1.0d0
+ ftprim(6)=0.0d0
+ ftbis(6)=0.0d0
+#endif
+ else if (rescale_mode.eq.2) then
+ quot=betaT/T0
+ quotl=1.0d0
+ do l=1,5
+ quotl1=quotl
+ quotl=quotl*quot
+ eplus=dexp(quotl)
+ eminus=dexp(-quotl)
+ logfac=1.0d0/dlog(eplus+eminus)
+ tanhT=(eplus-eminus)/(eplus+eminus)
+ fT(l)=1.12692801104297249644d0*logfac
+ ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
+ ftbis(l)=(l-1)*ftprim(l)/(quot*T0)-
+ & 2*l*quotl1/T0*logfac*
+ & (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2)
+ & +ftprim(l)*tanhT)
+ enddo
+#if defined(FUNCTH)
+ ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
+ & 320.0d0
+ ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
+ ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
+ & /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
+#elif defined(FUNCT)
+ fT(6)=betaT/T0
+ ftprim(6)=1.0d0/T0
+ ftbis(6)=0.0d0
+#else
+ fT(6)=1.0d0
+ ftprim(6)=0.0d0
+ ftbis(6)=0.0d0
+#endif
+ else if (rescale_mode.eq.0) then
+ do l=1,5
+ fT(l)=1.0d0
+ ftprim(l)=0.0d0
+ enddo
+ else
+ write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
+ & rescale_mode
+ call flush(iout)
+ return1
+ endif
+c write (iout,*) "ftprim",ftprim
+c write (iout,*) "ftbis",ftbis
+ betaT=1.0d0/(1.987D-3*betaT)
+#ifdef SPLITELE
+ etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
+ & +wvdwpp*evdw1
+ & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
+ & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
+ & +ft(2)*wturn3*eello_turn3
+ & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
+ & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
+ & +wbond*estr
+ eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees
+ & +ftprim(1)*wtor*etors+
+ & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
+ & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
+ & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
+ & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
+ & ftprim(1)*wsccor*esccor
+ ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+
+ & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
+ & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
+ & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
+ & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
+ & ftbis(1)*wsccor*esccor
+#else
+ etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
+ & +ft(1)*welec*(ees+evdw1)
+ & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
+ & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
+ & +ft(2)*wturn3*eello_turn3
+ & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
+ & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
+ & +wbond*estr
+ eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1)
+ & +ftprim(1)*wtor*etors+
+ & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
+ & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
+ & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
+ & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
+ & ftprim(1)*wsccor*esccor
+ ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
+ & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
+ & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
+ & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
+ & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
+ & ftprim(1)*wsccor*esccor
+#endif
+ weight=dexp(-betaT*(etot-potEmin)+entfac(t))
+#define DEBUG
+#ifdef DEBUG
+ write (iout,*) "iparm",iparm," t",t," betaT",betaT,
+ & " etot",etot," entfac",entfac(t),
+ & " weight",weight," ebis",ebis
+#endif
+#undef DEBUG
+ etot=etot-temper*eprim
+#ifdef MPI
+ sumW_p(k,iparm)=sumW_p(k,iparm)+weight
+ sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
+ sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
+ sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
+ do j=1,nQ+2
+ sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
+ sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
+ sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm)
+ & +etot*q(j,t)*weight
+ enddo
+#else
+ sumW(k,iparm)=sumW(k,iparm)+weight
+ sumE(k,iparm)=sumE(k,iparm)+etot*weight
+ sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
+ sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
+ do j=1,nQ+2
+ sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
+ sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
+ sumEQ(j,k,iparm)=sumEQ(j,k,iparm)
+ & +etot*q(j,t)*weight
+ enddo
+#endif
+ enddo
+ indE = aint(potE(t,iparm)-aint(potEmin))
+ if (indE.ge.0 .and. indE.le.maxinde) then
+ if (indE.gt.upindE_p) upindE_p=indE
+ histE_p(indE)=histE_p(indE)+dexp(-entfac(t))
+ endif
+#ifdef MPI
+ do ib=1,nT_h(iparm)
+ expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
+ hfin_p(ind,ib)=hfin_p(ind,ib)+
+ & dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
+ if (rmsrgymap) then
+ indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
+ indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
+ hrmsrgy_p(indrgy,indrms,ib)=
+ & hrmsrgy_p(indrgy,indrms,ib)+expfac
+ endif
+ enddo
+#else
+ do ib=1,nT_h(iparm)
+ expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
+ hfin(ind,ib)=hfin(ind,ib)+
+ & dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
+ if (rmsrgymap) then
+ indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
+ indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
+ hrmsrgy(indrgy,indrms,ib)=
+ & hrmsrgy(indrgy,indrms,ib)+expfac
+ endif
+ enddo
+#endif
+ enddo ! t
+ do ib=1,nT_h(iparm)
+ if (histout) call MPI_Reduce(hfin_p(0,ib),hfin(0,ib),nbin,
+ & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
+ if (rmsrgymap) then
+ call MPI_Reduce(hrmsrgy_p(0,0,ib),hrmsrgy(0,0,ib),
+ & (MaxBinRgy+1)*(nbin_rms+1),MPI_DOUBLE_PRECISION,MPI_SUM,Master,
+ & WHAM_COMM,IERROR)
+ endif
+ enddo
+ call MPI_Reduce(upindE_p,upindE,1,
+ & MPI_INTEGER,MPI_MAX,Master,WHAM_COMM,IERROR)
+ call MPI_Reduce(histE_p(0),histE(0),maxindE,
+ & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
+
+ if (me1.eq.master) then
+
+ if (histout) then
+
+ write (iout,'(6x,$)')
+ write (iout,'(f20.2,$)') (1.0d0/(1.987D-3*beta_h(ib,iparm)),
+ & ib=1,nT_h(iparm))
+ write (iout,*)
+
+ write (iout,'(/a)') 'Final histograms'
+ if (histfile) then
+ if (nslice.eq.1) then
+ if (separate_parset) then
+ write(licz3,"(bz,i3.3)") myparm
+ histname=prefix(:ilen(prefix))//'_par'//licz3//'.hist'
+ else
+ histname=prefix(:ilen(prefix))//'.hist'
+ endif
+ else
+ if (separate_parset) then
+ write(licz3,"(bz,i3.3)") myparm
+ histname=prefix(:ilen(prefix))//'_par'//licz3//
+ & '_slice_'//licz2//'.hist'
+ else
+ histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.hist'
+ endif
+ endif
+#if defined(AIX) || defined(PGI)
+ open (ihist,file=histname,position='append')
+#else
+ open (ihist,file=histname,access='append')
+#endif
+ endif
+
+ do t=0,tmax
+ liczba=t
+ sumH=0.0d0
+ do ib=1,nT_h(iparm)
+ sumH=sumH+hfin(t,ib)
+ enddo
+ if (sumH.gt.0.0d0) then
+ do j=1,nQ
+ jj = mod(liczba,nbin1)
+ liczba=liczba/nbin1
+ write (iout,'(f6.3,$)') dmin+(jj+0.5d0)*delta
+ if (histfile)
+ & write (ihist,'(f6.3,$)') dmin+(jj+0.5d0)*delta
+ enddo
+ do ib=1,nT_h(iparm)
+ write (iout,'(e20.10,$)') hfin(t,ib)
+ if (histfile) write (ihist,'(e20.10,$)') hfin(t,ib)
+ enddo
+ write (iout,'(i5)') iparm
+ if (histfile) write (ihist,'(i5)') iparm
+ endif
+ enddo
+
+ endif
+
+ if (entfile) then
+ if (nslice.eq.1) then
+ if (separate_parset) then
+ write(licz3,"(bz,i3.3)") myparm
+ histname=prefix(:ilen(prefix))//"_par"//licz3//'.ent'
+ else
+ histname=prefix(:ilen(prefix))//'.ent'
+ endif
+ else
+ if (separate_parset) then
+ write(licz3,"(bz,i3.3)") myparm
+ histname=prefix(:ilen(prefix))//'par_'//licz3//
+ & '_slice_'//licz2//'.ent'
+ else
+ histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.ent'
+ endif
+ endif
+#if defined(AIX) || defined(PGI)
+ open (ihist,file=histname,position='append')
+#else
+ open (ihist,file=histname,access='append')
+#endif
+ write (ihist,'(a)') "# Microcanonical entropy"
+ do i=0,upindE
+ write (ihist,'(f8.0,$)') dint(potEmin)+i
+ if (histE(i).gt.0.0e0) then
+ write (ihist,'(f15.5,$)') dlog(histE(i))
+ else
+ write (ihist,'(f15.5,$)') 0.0d0
+ endif
+ enddo
+ write (ihist,*)
+ close(ihist)
+ endif
+ write (iout,*) "Microcanonical entropy"
+ do i=0,upindE
+ write (iout,'(f8.0,$)') dint(potEmin)+i
+ if (histE(i).gt.0.0e0) then
+ write (iout,'(f15.5,$)') dlog(histE(i))
+ else
+ write (iout,'(f15.5,$)') 0.0d0
+ endif
+ write (iout,*)
+ enddo
+ if (rmsrgymap) then
+ if (nslice.eq.1) then
+ if (separate_parset) then
+ write(licz3,"(bz,i3.3)") myparm
+ histname=prefix(:ilen(prefix))//'_par'//licz3//'.rmsrgy'
+ else
+ histname=prefix(:ilen(prefix))//'.rmsrgy'
+ endif
+ else
+ if (separate_parset) then
+ write(licz3,"(bz,i3.3)") myparm
+ histname=prefix(:ilen(prefix))//'_par'//licz3//
+ & '_slice_'//licz2//'.rmsrgy'
+ else
+ histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.rmsrgy'
+ endif
+ endif
+#if defined(AIX) || defined(PGI)
+ open (ihist,file=histname,position='append')
+#else
+ open (ihist,file=histname,access='append')
+#endif
+ do i=0,nbin_rms
+ do j=0,nbin_rgy
+ write(ihist,'(2f8.2,$)')
+ & rgymin+deltrgy*j,rmsmin+deltrms*i
+ do ib=1,nT_h(iparm)
+ if (hrmsrgy(j,i,ib).gt.0.0d0) then
+ write(ihist,'(e14.5,$)')
+ & -dlog(hrmsrgy(j,i,ib))/beta_h(ib,iparm)
+ & +potEmin
+ else
+ write(ihist,'(e14.5,$)') 1.0d6
+ endif
+ enddo
+ write (ihist,'(i2)') iparm
+ enddo
+ enddo
+ close(ihist)
+ endif
+ endif
+ enddo ! iparm
+#ifdef MPI
+ call MPI_Reduce(hfin_ent_p(0),hfin_ent(0),nbin,
+ & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
+ call MPI_Reduce(sumW_p(0,1),sumW(0,1),(nGridT+1)*nParmSet,
+ & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
+ call MPI_Reduce(sumE_p(0,1),sumE(0,1),(nGridT+1)*nParmSet,
+ & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
+ call MPI_Reduce(sumEbis_p(0,1),sumEbis(0,1),(nGridT+1)*nParmSet,
+ & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
+ call MPI_Reduce(sumEsq_p(0,1),sumEsq(0,1),(nGridT+1)*nParmSet,
+ & MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
+ call MPI_Reduce(sumQ_p(1,0,1),sumQ(1,0,1),
+ & MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
+ & WHAM_COMM,IERROR)
+ call MPI_Reduce(sumQsq_p(1,0,1),sumQsq(1,0,1),
+ & MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
+ & WHAM_COMM,IERROR)
+ call MPI_Reduce(sumEQ_p(1,0,1),sumEQ(1,0,1),
+ & MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
+ & WHAM_COMM,IERROR)
+ if (me.eq.master) then
+#endif
+ write (iout,'(/a)') 'Thermal characteristics of folding'
+ if (nslice.eq.1) then
+ nazwa=prefix
+ else
+ nazwa=prefix(:ilen(prefix))//"_slice_"//licz2
+ endif
+ iln=ilen(nazwa)
+ if (nparmset.eq.1 .and. .not.separate_parset) then
+ nazwa=nazwa(:iln)//".thermal"
+ else if (nparmset.eq.1 .and. separate_parset) then
+ write(licz3,"(bz,i3.3)") myparm
+ nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
+ endif
+ do iparm=1,nParmSet
+ if (nparmset.gt.1) then
+ write(licz3,"(bz,i3.3)") iparm
+ nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
+ endif
+ open(34,file=nazwa)
+ if (separate_parset) then
+ write (iout,'(a,i3)') "Parameter set",myparm
+ else
+ write (iout,'(a,i3)') "Parameter set",iparm
+ endif
+ do i=0,NGridT
+ sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
+ sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/
+ & sumW(i,iparm)
+ sumEsq(i,iparm)=(sumEsq(i,iparm)/sumW(i,iparm)
+ & -sumE(i,iparm)**2)/(1.987D-3*(startGridT+i*delta_T)**2)
+ do j=1,nQ+2
+ sumQ(j,i,iparm)=sumQ(j,i,iparm)/sumW(i,iparm)
+ sumQsq(j,i,iparm)=sumQsq(j,i,iparm)/sumW(i,iparm)
+ & -sumQ(j,i,iparm)**2
+ sumEQ(j,i,iparm)=sumEQ(j,i,iparm)/sumW(i,iparm)
+ & -sumQ(j,i,iparm)*sumE(i,iparm)
+ enddo
+ sumW(i,iparm)=-dlog(sumW(i,iparm))*(1.987D-3*
+ & (startGridT+i*delta_T))+potEmin
+ write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
+ & sumW(i,iparm),sumE(i,iparm)
+ write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
+ write (iout,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
+ & (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
+ write (iout,*)
+ write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
+ & sumW(i,iparm),sumE(i,iparm)
+ write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
+ write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
+ & (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
+ write (34,*)
+ enddo
+ close(34)
+ enddo
+ if (histout) then
+ do t=0,tmax
+ if (hfin_ent(t).gt.0.0d0) then
+ liczba=t
+ jj = mod(liczba,nbin1)
+ write (iout,'(f6.3,e20.10," ent")') dmin+(jj+0.5d0)*delta,
+ & hfin_ent(t)
+ if (histfile) write (ihist,'(f6.3,e20.10," ent")')
+ & dmin+(jj+0.5d0)*delta,
+ & hfin_ent(t)
+ endif
+ enddo
+ if (histfile) close(ihist)
+ endif
+
+#ifdef ZSCORE
+! Write data for zscore
+ if (nslice.eq.1) then
+ zscname=prefix(:ilen(prefix))//".zsc"
+ else
+ zscname=prefix(:ilen(prefix))//"_slice_"//licz2//".zsc"
+ endif
+#if defined(AIX) || defined(PGI)
+ open (izsc,file=prefix(:ilen(prefix))//'.zsc',position='append')
+#else
+ open (izsc,file=prefix(:ilen(prefix))//'.zsc',access='append')
+#endif
+ write (izsc,'("NQ=",i1," NPARM=",i1)') nQ,nParmSet
+ do iparm=1,nParmSet
+ write (izsc,'("NT=",i1)') nT_h(iparm)
+ do ib=1,nT_h(iparm)
+ write (izsc,'("TEMP=",f6.1," NR=",i2," SNK=",$)')
+ & 1.0d0/(beta_h(ib,iparm)*1.987D-3),nR(ib,iparm)
+ jj = min0(nR(ib,iparm),7)
+ write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=1,jj)
+ write (izsc,'(a1,$)') (" ",i=22+8*jj+1,79)
+ write (izsc,'("&")')
+ if (nR(ib,iparm).gt.7) then
+ do ii=8,nR(ib,iparm),9
+ jj = min0(nR(ib,iparm),ii+8)
+ write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=ii,jj)
+ write (izsc,'(a1,$') (" ",i=(jj-ii+1)*8+1,79)
+ write (izsc,'("&")')
+ enddo
+ endif
+ write (izsc,'("FI=",$)')
+ jj=min0(nR(ib,iparm),7)
+ write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=1,jj)
+ write (izsc,'(a1,$)') (" ",i=3+10*jj+1,79)
+ write (izsc,'("&")')
+ if (nR(ib,iparm).gt.7) then
+ do ii=8,nR(ib,iparm),9
+ jj = min0(nR(ib,iparm),ii+8)
+ write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=ii,jj)
+ if (jj.eq.nR(ib,iparm)) then
+ write (izsc,*)
+ else
+ write (izsc,'(a1,$)') (" ",i=10*(jj-ii+1)+1,79)
+ write (izsc,'(t80,"&")')
+ endif
+ enddo
+ endif
+ do i=1,nR(ib,iparm)
+ write (izsc,'("KH=",$)')
+ write (izsc,'(f7.2,$)') (Kh(j,i,ib,iparm),j=1,nQ)
+ write (izsc,'(" Q0=",$)')
+ write (izsc,'(f7.5,$)') (q0(j,i,ib,iparm),j=1,nQ)
+ write (izsc,*)
+ enddo
+ enddo
+ enddo
+ close(izsc)
+#endif
+#ifdef MPI
+ endif
+#endif
+
+ return
+
+ end
projects / unres.git / commitdiff