--- /dev/null
+ module io_wham
+
+ use io_units
+ use io_base
+ use wham_data
+#ifndef CLUSTER
+ use w_compar_data
+#endif
+! use geometry_data
+! use geometry
+ implicit none
+!-----------------------------------------------------------------------------
+!
+!
+!-----------------------------------------------------------------------------
+ contains
+!-----------------------------------------------------------------------------
+! openunits.F
+!-----------------------------------------------------------------------------
+#ifndef CLUSTER
+ subroutine openunits
+#ifdef WIN
+ use dfport
+#endif
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'DIMENSIONS.ZSCOPT'
+#ifdef MPI
+ use MPI_data
+ include 'mpif.h'
+! include 'COMMON.MPI'
+! integer :: MyRank
+ character(len=3) :: liczba
+#endif
+! include 'COMMON.IOUNITS'
+ integer :: lenpre,lenpot !,ilen
+!el external ilen
+
+#ifdef MPI
+ MyRank=Me
+#endif
+ call mygetenv('PREFIX',prefix)
+ call mygetenv('SCRATCHDIR',scratchdir)
+ call mygetenv('POT',pot)
+ lenpre=ilen(prefix)
+ lenpot=ilen(pot)
+ call mygetenv('POT',pot)
+ entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
+! Get the names and open the input files
+ open (1,file=prefix(:ilen(prefix))//'.inp',status='old')
+! Get parameter filenames and open the parameter files.
+ call mygetenv('BONDPAR',bondname)
+ open (ibond,file=bondname,status='old')
+ call mygetenv('THETPAR',thetname)
+ open (ithep,file=thetname,status='old')
+ call mygetenv('ROTPAR',rotname)
+ open (irotam,file=rotname,status='old')
+ call mygetenv('TORPAR',torname)
+ open (itorp,file=torname,status='old')
+ call mygetenv('TORDPAR',tordname)
+ open (itordp,file=tordname,status='old')
+ call mygetenv('FOURIER',fouriername)
+ open (ifourier,file=fouriername,status='old')
+ call mygetenv('SCCORPAR',sccorname)
+ open (isccor,file=sccorname,status='old')
+ call mygetenv('ELEPAR',elename)
+ open (ielep,file=elename,status='old')
+ call mygetenv('SIDEPAR',sidename)
+ open (isidep,file=sidename,status='old')
+ call mygetenv('SIDEP',sidepname)
+ open (isidep1,file=sidepname,status="old")
+#ifndef OLDSCP
+!
+! 8/9/01 In the newest version SCp interaction constants are read from a file
+! Use -DOLDSCP to use hard-coded constants instead.
+!
+ call mygetenv('SCPPAR',scpname)
+ open (iscpp,file=scpname,status='old')
+#endif
+#ifdef MPL
+ if (MyID.eq.BossID) then
+ MyRank = MyID/fgProcs
+#endif
+#ifdef MPI
+ print *,'OpenUnits: processor',MyRank
+ call numstr(MyRank,liczba)
+ outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//liczba
+#else
+ outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
+#endif
+ open(iout,file=outname,status='unknown')
+ write (iout,'(80(1h-))')
+ write (iout,'(30x,a)') "FILE ASSIGNMENT"
+ write (iout,'(80(1h-))')
+ write (iout,*) "Input file : ",&
+ prefix(:ilen(prefix))//'.inp'
+ write (iout,*) "Output file : ",&
+ outname(:ilen(outname))
+ write (iout,*)
+ write (iout,*) "Sidechain potential file : ",&
+ sidename(:ilen(sidename))
+#ifndef OLDSCP
+ write (iout,*) "SCp potential file : ",&
+ scpname(:ilen(scpname))
+#endif
+ write (iout,*) "Electrostatic potential file : ",&
+ elename(:ilen(elename))
+ write (iout,*) "Cumulant coefficient file : ",&
+ fouriername(:ilen(fouriername))
+ write (iout,*) "Torsional parameter file : ",&
+ torname(:ilen(torname))
+ write (iout,*) "Double torsional parameter file : ",&
+ tordname(:ilen(tordname))
+ write (iout,*) "Backbone-rotamer parameter file : ",&
+ sccorname(:ilen(sccorname))
+ write (iout,*) "Bond & inertia constant file : ",&
+ bondname(:ilen(bondname))
+ write (iout,*) "Bending parameter file : ",&
+ thetname(:ilen(thetname))
+ write (iout,*) "Rotamer parameter file : ",&
+ rotname(:ilen(rotname))
+ write (iout,'(80(1h-))')
+ write (iout,*)
+ return
+ end subroutine openunits
+!-----------------------------------------------------------------------------
+! molread_zs.F
+!-----------------------------------------------------------------------------
+ subroutine molread(*)
+!
+! Read molecular data.
+!
+ use energy_data
+ use geometry_data, only:nres,deg2rad,c,dc
+ use control_data, only:iscode
+ use control, only:rescode,setup_var,init_int_table
+ use geometry, only:alloc_geo_arrays
+ use energy, only:alloc_ener_arrays
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'DIMENSIONS.ZSCOPT'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.GEO'
+! include 'COMMON.VAR'
+! include 'COMMON.INTERACT'
+! include 'COMMON.LOCAL'
+! include 'COMMON.NAMES'
+! include 'COMMON.CHAIN'
+! include 'COMMON.FFIELD'
+! include 'COMMON.SBRIDGE'
+! include 'COMMON.TORCNSTR'
+! include 'COMMON.CONTROL'
+ character(len=4),dimension(:),allocatable :: sequence !(nres)
+!el integer :: rescode
+!el real(kind=8) :: x(maxvar)
+ character(len=320) :: controlcard !,ucase
+ integer,dimension(nres) :: itype_pdb !(maxres)
+ integer :: i,j,i1,i2,it1,it2
+ real(kind=8) :: scalscp
+!el logical :: seq_comp
+ call card_concat(controlcard,.true.)
+ call reada(controlcard,'SCAL14',scal14,0.4d0)
+ call reada(controlcard,'SCALSCP',scalscp,1.0d0)
+ call reada(controlcard,'CUTOFF',cutoff_corr,7.0d0)
+ call reada(controlcard,'TEMP0',temp0,300.0d0) !el
+ call reada(controlcard,'DELT_CORR',delt_corr,0.5d0)
+ r0_corr=cutoff_corr-delt_corr
+ call readi(controlcard,"NRES",nres,0)
+ allocate(sequence(nres+1))
+!el znamy juz ilosc reszt wiec mozna zaalokowac tablice do liczenia enerii
+ call alloc_geo_arrays
+ call alloc_ener_arrays
+! alokacja dodatkowych tablic, ktore w unresie byly alokowanie w locie
+!----------------------------
+ allocate(c(3,2*nres+2))
+ allocate(dc(3,0:2*nres+2))
+ allocate(itype(nres+2))
+ allocate(itel(nres+2))
+!
+! Zero out tableis.
+ do i=1,2*nres+2
+ do j=1,3
+ c(j,i)=0.0D0
+ dc(j,i)=0.0D0
+ enddo
+ enddo
+ do i=1,nres+2
+ itype(i)=0
+ itel(i)=0
+ enddo
+!--------------------------
+!
+ iscode=index(controlcard,"ONE_LETTER")
+ if (nres.le.0) then
+ write (iout,*) "Error: no residues in molecule"
+ return 1
+ endif
+ if (nres.gt.maxres) then
+ write (iout,*) "Error: too many residues",nres,maxres
+ endif
+ write(iout,*) 'nres=',nres
+! Read sequence of the protein
+ if (iscode.gt.0) then
+ read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
+ else
+ read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
+ endif
+! Convert sequence to numeric code
+ do i=1,nres
+ itype(i)=rescode(i,sequence(i),iscode)
+ enddo
+ write (iout,*) "Numeric code:"
+ write (iout,'(20i4)') (itype(i),i=1,nres)
+ do i=1,nres-1
+#ifdef PROCOR
+ if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
+#else
+ if (itype(i).eq.ntyp1) then
+#endif
+ itel(i)=0
+#ifdef PROCOR
+ else if (iabs(itype(i+1)).ne.20) then
+#else
+ else if (iabs(itype(i)).ne.20) then
+#endif
+ itel(i)=1
+ else
+ itel(i)=2
+ endif
+ enddo
+ write (iout,*) "ITEL"
+ do i=1,nres-1
+ write (iout,*) i,itype(i),itel(i)
+ enddo
+ 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
+ if (itype(1).eq.ntyp1) nnt=2
+ if (itype(nres).eq.ntyp1) nct=nct-1
+ write(iout,*) 'NNT=',NNT,' NCT=',NCT
+ call setup_var
+ call init_int_table
+ if (ns.gt.0) then
+ write (iout,'(/a,i3,a)') 'The chain contains',ns,&
+ ' disulfide-bridging cysteines.'
+ write (iout,'(20i4)') (iss(i),i=1,ns)
+ write (iout,'(/a/)') 'Pre-formed links are:'
+ do i=1,nss
+ i1=ihpb(i)-nres
+ i2=jhpb(i)-nres
+ it1=itype(i1)
+ it2=itype(i2)
+ write (iout,'(2a,i3,3a,i3,a,3f10.3)') &
+ restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',&
+ dhpb(i),ebr,forcon(i)
+ enddo
+ endif
+ write (iout,'(a)')
+ return
+ end subroutine molread
+!-----------------------------------------------------------------------------
+! parmread.F
+!-----------------------------------------------------------------------------
+ subroutine parmread(iparm,*)
+#else
+ subroutine parmread
+#endif
+!
+! Read the parameters of the probability distributions of the virtual-bond
+! valence angles and the side chains and energy parameters.
+!
+ use wham_data
+
+ use geometry_data
+ use energy_data
+ use control_data, only: maxtor,maxterm,maxlor,maxterm_sccor,&
+ maxtermd_1,maxtermd_2,maxthetyp,maxthetyp1
+ use MD_data
+!el use MPI_data
+!el use map_data
+ use io_config, only: printmat
+ use control, only: getenv_loc
+
+#ifdef MPI
+ use MPI_data
+ include "mpif.h"
+ integer :: IERROR
+#endif
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'DIMENSIONS.ZSCOPT'
+! include 'DIMENSIONS.FREE'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CHAIN'
+! include 'COMMON.INTERACT'
+! include 'COMMON.GEO'
+! include 'COMMON.LOCAL'
+! include 'COMMON.TORSION'
+! include 'COMMON.FFIELD'
+! include 'COMMON.NAMES'
+! include 'COMMON.SBRIDGE'
+! include 'COMMON.WEIGHTS'
+! include 'COMMON.ENEPS'
+! include 'COMMON.SCCOR'
+! include 'COMMON.SCROT'
+! include 'COMMON.FREE'
+ character(len=1) :: t1,t2,t3
+ character(len=1) :: onelett(4) = (/"G","A","P","D"/)
+ character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
+ logical :: lprint
+ real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob)
+ character(len=800) :: controlcard
+ character(len=256) :: bondname_t,thetname_t,rotname_t,torname_t,&
+ tordname_t,fouriername_t,elename_t,sidename_t,scpname_t,&
+ sccorname_t
+!el integer ilen
+!el external ilen
+ character(len=16) :: key
+ integer :: iparm
+!el real(kind=8) :: ip,mp
+ real(kind=8) :: dwa16,akl,si,rri,epsij,rrij,sigeps,sigt1sq,&
+ sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm
+ real(kind=8) :: v0ij,v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,rjunk,&
+ res1
+ integer :: i,j,ichir1,ichir2,k,l,m,kk,ii,mm,junk,lll,ll,llll,n
+ integer :: nlobi,iblock,maxinter,iscprol
+!
+! Body
+!
+! Set LPRINT=.TRUE. for debugging
+ dwa16=2.0d0**(1.0d0/6.0d0)
+ lprint=.false.
+ itypro=20
+! Assign virtual-bond length
+ vbl=3.8D0
+ vblinv=1.0D0/vbl
+ vblinv2=vblinv*vblinv
+#ifndef CLUSTER
+ call card_concat(controlcard,.true.)
+ wname(4)="WCORRH"
+!el
+allocate(ww(max_eneW))
+ do i=1,n_eneW
+ key = wname(i)(:ilen(wname(i)))
+ call reada(controlcard,key(:ilen(key)),ww(i),1.0d0)
+ enddo
+
+ write (iout,*) "iparm",iparm," myparm",myparm
+! If reading not own parameters, skip assignment
+
+ if (iparm.eq.myparm .or. .not.separate_parset) then
+
+!
+! Setup weights for UNRES
+!
+ wsc=ww(1)
+ wscp=ww(2)
+ welec=ww(3)
+ wcorr=ww(4)
+ wcorr5=ww(5)
+ wcorr6=ww(6)
+ wel_loc=ww(7)
+ wturn3=ww(8)
+ wturn4=ww(9)
+ wturn6=ww(10)
+ wang=ww(11)
+ wscloc=ww(12)
+ wtor=ww(13)
+ wtor_d=ww(14)
+ wvdwpp=ww(16)
+ wbond=ww(18)
+ wsccor=ww(19)
+
+ endif
+!
+!el------
+ allocate(weights(n_ene))
+ weights(1)=wsc
+ weights(2)=wscp
+ weights(3)=welec
+ weights(4)=wcorr
+ weights(5)=wcorr5
+ weights(6)=wcorr6
+ weights(7)=wel_loc
+ weights(8)=wturn3
+ weights(9)=wturn4
+ weights(10)=wturn6
+ weights(11)=wang
+ weights(12)=wscloc
+ weights(13)=wtor
+ weights(14)=wtor_d
+ weights(15)=0 !wstrain !
+ weights(16)=0 !wvdwpp !
+ weights(17)=wbond
+ weights(18)=0 !scal14 !
+ weights(21)=wsccor
+! el--------
+ call card_concat(controlcard,.false.)
+
+! Return if not own parameters
+
+ if (iparm.ne.myparm .and. separate_parset) return
+
+ call reads(controlcard,"BONDPAR",bondname_t,bondname)
+ open (ibond,file=bondname_t,status='old')
+ rewind(ibond)
+ call reads(controlcard,"THETPAR",thetname_t,thetname)
+ open (ithep,file=thetname_t,status='old')
+ rewind(ithep)
+ call reads(controlcard,"ROTPAR",rotname_t,rotname)
+ open (irotam,file=rotname_t,status='old')
+ rewind(irotam)
+ call reads(controlcard,"TORPAR",torname_t,torname)
+ open (itorp,file=torname_t,status='old')
+ rewind(itorp)
+ call reads(controlcard,"TORDPAR",tordname_t,tordname)
+ open (itordp,file=tordname_t,status='old')
+ rewind(itordp)
+ call reads(controlcard,"SCCORPAR",sccorname_t,sccorname)
+ open (isccor,file=sccorname_t,status='old')
+ rewind(isccor)
+ call reads(controlcard,"FOURIER",fouriername_t,fouriername)
+ open (ifourier,file=fouriername_t,status='old')
+ rewind(ifourier)
+ call reads(controlcard,"ELEPAR",elename_t,elename)
+ open (ielep,file=elename_t,status='old')
+ rewind(ielep)
+ call reads(controlcard,"SIDEPAR",sidename_t,sidename)
+ open (isidep,file=sidename_t,status='old')
+ rewind(isidep)
+ call reads(controlcard,"SCPPAR",scpname_t,scpname)
+ open (iscpp,file=scpname_t,status='old')
+ rewind(iscpp)
+ write (iout,*) "Parameter set:",iparm
+ write (iout,*) "Energy-term weights:"
+ do i=1,n_eneW
+ write (iout,'(a16,f10.5)') wname(i),ww(i)
+ enddo
+ write (iout,*) "Sidechain potential file : ",&
+ sidename_t(:ilen(sidename_t))
+#ifndef OLDSCP
+ write (iout,*) "SCp potential file : ",&
+ scpname_t(:ilen(scpname_t))
+#endif
+ write (iout,*) "Electrostatic potential file : ",&
+ elename_t(:ilen(elename_t))
+ write (iout,*) "Cumulant coefficient file : ",&
+ fouriername_t(:ilen(fouriername_t))
+ write (iout,*) "Torsional parameter file : ",&
+ torname_t(:ilen(torname_t))
+ write (iout,*) "Double torsional parameter file : ",&
+ tordname_t(:ilen(tordname_t))
+ write (iout,*) "Backbone-rotamer parameter file : ",&
+ sccorname(:ilen(sccorname))
+ write (iout,*) "Bond & inertia constant file : ",&
+ bondname_t(:ilen(bondname_t))
+ write (iout,*) "Bending parameter file : ",&
+ thetname_t(:ilen(thetname_t))
+ write (iout,*) "Rotamer parameter file : ",&
+ rotname_t(:ilen(rotname_t))
+#endif
+!
+! Read the virtual-bond parameters, masses, and moments of inertia
+! and Stokes' radii of the peptide group and side chains
+!
+ allocate(dsc(ntyp1)) !(ntyp1)
+ allocate(dsc_inv(ntyp1)) !(ntyp1)
+ allocate(nbondterm(ntyp)) !(ntyp)
+ allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
+ allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
+!el allocate(msc(ntyp+1)) !(ntyp+1)
+!el allocate(isc(ntyp+1)) !(ntyp+1)
+!el allocate(restok(ntyp+1)) !(ntyp+1)
+ allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
+
+#ifdef CRYST_BOND
+ read (ibond,*) vbldp0,akp
+ do i=1,ntyp
+ nbondterm(i)=1
+ read (ibond,*) vbldsc0(1,i),aksc(1,i)
+ dsc(i) = vbldsc0(1,i)
+ if (i.eq.10) then
+ dsc_inv(i)=0.0D0
+ else
+ dsc_inv(i)=1.0D0/dsc(i)
+ endif
+ enddo
+#else
+ read (ibond,*) ijunk,vbldp0,akp,rjunk
+ do i=1,ntyp
+ read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
+ j=1,nbondterm(i))
+ dsc(i) = vbldsc0(1,i)
+ if (i.eq.10) then
+ dsc_inv(i)=0.0D0
+ else
+ dsc_inv(i)=1.0D0/dsc(i)
+ endif
+ enddo
+#endif
+ if (lprint) then
+ write(iout,'(/a/)')"Force constants virtual bonds:"
+ write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',&
+ 'inertia','Pstok'
+ write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
+ do i=1,ntyp
+ write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),&
+ vbldsc0(1,i),aksc(1,i),abond0(1,i)
+ do j=2,nbondterm(i)
+ write (iout,'(13x,3f10.5)') &
+ vbldsc0(j,i),aksc(j,i),abond0(j,i)
+ enddo
+ enddo
+ endif
+!----------------------------------------------------
+ allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
+ allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
+ allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
+ allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
+ allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
+ allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
+ do i=-ntyp,ntyp
+ a0thet(i)=0.0D0
+ do j=1,2
+ do ichir1=-1,1
+ do ichir2=-1,1
+ athet(j,i,ichir1,ichir2)=0.0D0
+ bthet(j,i,ichir1,ichir2)=0.0D0
+ enddo
+ enddo
+ enddo
+ do j=0,3
+ polthet(j,i)=0.0D0
+ enddo
+ do j=1,3
+ gthet(j,i)=0.0D0
+ enddo
+ theta0(i)=0.0D0
+ sig0(i)=0.0D0
+ sigc0(i)=0.0D0
+ enddo
+!elwrite(iout,*) "parmread kontrol"
+
+#ifdef CRYST_THETA
+!
+! Read the parameters of the probability distribution/energy expression
+! of the virtual-bond valence angles theta
+!
+ do i=1,ntyp
+ read (ithep,*) a0thet(i),(athet(j,i,1,1),j=1,2),&
+ (bthet(j,i,1,1),j=1,2)
+ read (ithep,*) (polthet(j,i),j=0,3)
+!elwrite(iout,*) "parmread kontrol in cryst_theta"
+ read (ithep,*) (gthet(j,i),j=1,3)
+!elwrite(iout,*) "parmread kontrol in cryst_theta"
+ read (ithep,*) theta0(i),sig0(i),sigc0(i)
+ sigc0(i)=sigc0(i)**2
+!elwrite(iout,*) "parmread kontrol in cryst_theta"
+ enddo
+!elwrite(iout,*) "parmread kontrol in cryst_theta"
+ do i=1,ntyp
+ athet(1,i,1,-1)=athet(1,i,1,1)
+ athet(2,i,1,-1)=athet(2,i,1,1)
+ bthet(1,i,1,-1)=-bthet(1,i,1,1)
+ bthet(2,i,1,-1)=-bthet(2,i,1,1)
+ athet(1,i,-1,1)=-athet(1,i,1,1)
+ athet(2,i,-1,1)=-athet(2,i,1,1)
+ bthet(1,i,-1,1)=bthet(1,i,1,1)
+ bthet(2,i,-1,1)=bthet(2,i,1,1)
+ enddo
+!elwrite(iout,*) "parmread kontrol in cryst_theta"
+ do i=-ntyp,-1
+ a0thet(i)=a0thet(-i)
+ athet(1,i,-1,-1)=athet(1,-i,1,1)
+ athet(2,i,-1,-1)=-athet(2,-i,1,1)
+ bthet(1,i,-1,-1)=bthet(1,-i,1,1)
+ bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
+ athet(1,i,-1,1)=athet(1,-i,1,1)
+ athet(2,i,-1,1)=-athet(2,-i,1,1)
+ bthet(1,i,-1,1)=-bthet(1,-i,1,1)
+ bthet(2,i,-1,1)=bthet(2,-i,1,1)
+ athet(1,i,1,-1)=-athet(1,-i,1,1)
+ athet(2,i,1,-1)=athet(2,-i,1,1)
+ bthet(1,i,1,-1)=bthet(1,-i,1,1)
+ bthet(2,i,1,-1)=-bthet(2,-i,1,1)
+ theta0(i)=theta0(-i)
+ sig0(i)=sig0(-i)
+ sigc0(i)=sigc0(-i)
+ do j=0,3
+ polthet(j,i)=polthet(j,-i)
+ enddo
+ do j=1,3
+ gthet(j,i)=gthet(j,-i)
+ enddo
+ enddo
+!elwrite(iout,*) "parmread kontrol in cryst_theta"
+ close (ithep)
+!elwrite(iout,*) "parmread kontrol in cryst_theta"
+ if (lprint) then
+! write (iout,'(a)')
+! & 'Parameters of the virtual-bond valence angles:'
+! write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
+! & ' ATHETA0 ',' A1 ',' A2 ',
+! & ' B1 ',' B2 '
+! do i=1,ntyp
+! write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
+! & a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
+! enddo
+! write (iout,'(/a/9x,5a/79(1h-))')
+! & 'Parameters of the expression for sigma(theta_c):',
+! & ' ALPH0 ',' ALPH1 ',' ALPH2 ',
+! & ' ALPH3 ',' SIGMA0C '
+! do i=1,ntyp
+! write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
+! & (polthet(j,i),j=0,3),sigc0(i)
+! enddo
+! write (iout,'(/a/9x,5a/79(1h-))')
+! & 'Parameters of the second gaussian:',
+! & ' THETA0 ',' SIGMA0 ',' G1 ',
+! & ' G2 ',' G3 '
+! do i=1,ntyp
+! write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
+! & sig0(i),(gthet(j,i),j=1,3)
+! enddo
+ write (iout,'(a)') &
+ 'Parameters of the virtual-bond valence angles:'
+ write (iout,'(/a/9x,5a/79(1h-))') &
+ 'Coefficients of expansion',&
+ ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
+ ' b1*10^1 ',' b2*10^1 '
+ do i=1,ntyp
+ write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),&
+ a0thet(i),(100*athet(j,i,1,1),j=1,2),&
+ (10*bthet(j,i,1,1),j=1,2)
+ enddo
+ write (iout,'(/a/9x,5a/79(1h-))') &
+ 'Parameters of the expression for sigma(theta_c):',&
+ ' alpha0 ',' alph1 ',' alph2 ',&
+ ' alhp3 ',' sigma0c '
+ do i=1,ntyp
+ write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),&
+ (polthet(j,i),j=0,3),sigc0(i)
+ enddo
+ write (iout,'(/a/9x,5a/79(1h-))') &
+ 'Parameters of the second gaussian:',&
+ ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
+ ' G2 ',' G3*10^1 '
+ do i=1,ntyp
+ write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),&
+ 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
+ enddo
+ endif
+#else
+!
+! Read the parameters of Utheta determined from ab initio surfaces
+! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
+!
+! write (iout,*) "tu dochodze"
+ read (ithep,*) nthetyp,ntheterm,ntheterm2,&
+ ntheterm3,nsingle,ndouble
+ nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
+
+!----------------------------------------------------
+ allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
+ allocate(aa0thet(-maxthetyp1:maxthetyp1,&
+ -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
+!(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
+ allocate(aathet(ntheterm,-maxthetyp1:maxthetyp1,&
+ -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
+!(maxtheterm,-maxthetyp1:maxthetyp1,&
+! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
+ allocate(bbthet(nsingle,ntheterm2,-maxthetyp1:maxthetyp1,&
+ -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
+ allocate(ccthet(nsingle,ntheterm2,-maxthetyp1:maxthetyp1,&
+ -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
+ allocate(ddthet(nsingle,ntheterm2,-maxthetyp1:maxthetyp1,&
+ -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
+ allocate(eethet(nsingle,ntheterm2,-maxthetyp1:maxthetyp1,&
+ -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
+!(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
+! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
+ allocate(ffthet(ndouble,ndouble,ntheterm3,-maxthetyp1:maxthetyp1,&
+ -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
+ allocate(ggthet(ndouble,ndouble,ntheterm3,-maxthetyp1:maxthetyp1,&
+ -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
+!(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
+! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
+
+
+ read (ithep,*) (ithetyp(i),i=1,ntyp1)
+ do i=-ntyp1,-1
+ ithetyp(i)=-ithetyp(-i)
+ enddo
+! write (iout,*) "tu dochodze"
+ do iblock=1,2
+ do i=-maxthetyp,maxthetyp
+ do j=-maxthetyp,maxthetyp
+ do k=-maxthetyp,maxthetyp
+ aa0thet(i,j,k,iblock)=0.0d0
+ do l=1,ntheterm
+ aathet(l,i,j,k,iblock)=0.0d0
+ enddo
+ do l=1,ntheterm2
+ do m=1,nsingle
+ bbthet(m,l,i,j,k,iblock)=0.0d0
+ ccthet(m,l,i,j,k,iblock)=0.0d0
+ ddthet(m,l,i,j,k,iblock)=0.0d0
+ eethet(m,l,i,j,k,iblock)=0.0d0
+ enddo
+ enddo
+ do l=1,ntheterm3
+ do m=1,ndouble
+ do mm=1,ndouble
+ ffthet(mm,m,l,i,j,k,iblock)=0.0d0
+ ggthet(mm,m,l,i,j,k,iblock)=0.0d0
+ enddo
+ enddo
+ enddo
+ enddo
+ enddo
+ enddo
+ enddo
+ do iblock=1,2
+ do i=0,nthetyp
+ do j=-nthetyp,nthetyp
+ do k=-nthetyp,nthetyp
+ read (ithep,'(6a)') res1
+ read (ithep,*) aa0thet(i,j,k,iblock)
+ read (ithep,*)(aathet(l,i,j,k,iblock),l=1,ntheterm)
+ read (ithep,*) &
+ ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
+ (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
+ (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
+ (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
+ ll=1,ntheterm2)
+ read (ithep,*) &
+ (((ffthet(llll,lll,ll,i,j,k,iblock),&
+ ffthet(lll,llll,ll,i,j,k,iblock),&
+ ggthet(llll,lll,ll,i,j,k,iblock),&
+ ggthet(lll,llll,ll,i,j,k,iblock),&
+ llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
+ enddo
+ enddo
+ enddo
+!
+! For dummy ends assign glycine-type coefficients of theta-only terms; the
+! coefficients of theta-and-gamma-dependent terms are zero.
+!
+ do i=1,nthetyp
+ do j=1,nthetyp
+ do l=1,ntheterm
+ aathet(l,i,j,nthetyp+1,iblock)=0.0d0
+ aathet(l,nthetyp+1,i,j,iblock)=0.0d0
+ enddo
+ aa0thet(i,j,nthetyp+1,iblock)=0.0d0
+ aa0thet(nthetyp+1,i,j,iblock)=0.0d0
+ enddo
+ do l=1,ntheterm
+ aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
+ enddo
+ aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
+ enddo
+ enddo
+! Substitution for D aminoacids from symmetry.
+ do iblock=1,2
+ do i=-nthetyp,0
+ do j=-nthetyp,nthetyp
+ do k=-nthetyp,nthetyp
+ aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
+ do l=1,ntheterm
+ aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
+ enddo
+ do ll=1,ntheterm2
+ do lll=1,nsingle
+ bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
+ ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
+ ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
+ eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
+ enddo
+ enddo
+ do ll=1,ntheterm3
+ do lll=2,ndouble
+ do llll=1,lll-1
+ ffthet(llll,lll,ll,i,j,k,iblock)= &
+ ffthet(llll,lll,ll,-i,-j,-k,iblock)
+ ffthet(lll,llll,ll,i,j,k,iblock)= &
+ ffthet(lll,llll,ll,-i,-j,-k,iblock)
+ ggthet(llll,lll,ll,i,j,k,iblock)= &
+ -ggthet(llll,lll,ll,-i,-j,-k,iblock)
+ ggthet(lll,llll,ll,i,j,k,iblock)= &
+ -ggthet(lll,llll,ll,-i,-j,-k,iblock)
+ enddo !ll
+ enddo !lll
+ enddo !llll
+ enddo !k
+ enddo !j
+ enddo !i
+ enddo !iblock
+
+!
+! Control printout of the coefficients of virtual-bond-angle potentials
+!
+do iblock=1,2
+ if (lprint) then
+ write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
+ do i=1,nthetyp+1
+ do j=1,nthetyp+1
+ do k=1,nthetyp+1
+ write (iout,'(//4a)') &
+ 'Type ',onelett(i),onelett(j),onelett(k)
+ write (iout,'(//a,10x,a)') " l","a[l]"
+ write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
+ write (iout,'(i2,1pe15.5)') &
+ (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
+ do l=1,ntheterm2
+ write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))') &
+ "b",l,"c",l,"d",l,"e",l
+ do m=1,nsingle
+ write (iout,'(i2,4(1pe15.5))') m,&
+ bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
+ ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
+ enddo
+ enddo
+ do l=1,ntheterm3
+ write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))') &
+ "f+",l,"f-",l,"g+",l,"g-",l
+ do m=2,ndouble
+ do n=1,m-1
+ write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
+ ffthet(n,m,l,i,j,k,iblock),&
+ ffthet(m,n,l,i,j,k,iblock),&
+ ggthet(n,m,l,i,j,k,iblock),&
+ ggthet(m,n,l,i,j,k,iblock)
+ enddo
+ enddo
+ enddo
+ enddo
+ enddo
+ enddo
+ call flush(iout)
+ endif
+enddo
+#endif
+!-------------------------------------------
+ allocate(nlob(ntyp1)) !(ntyp1)
+ allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
+ allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
+ allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
+
+ do i=1,ntyp
+ do j=1,maxlob
+ bsc(j,i)=0.0D0
+ nlob(i)=0
+ enddo
+ enddo
+ nlob(ntyp1)=0
+ dsc(ntyp1)=0.0D0
+
+ do i=-ntyp,ntyp
+ do j=1,maxlob
+ do k=1,3
+ censc(k,j,i)=0.0D0
+ enddo
+ do k=1,3
+ do l=1,3
+ gaussc(l,k,j,i)=0.0D0
+ enddo
+ enddo
+ enddo
+ enddo
+
+#ifdef CRYST_SC
+!
+! Read the parameters of the probability distribution/energy expression
+! of the side chains.
+!
+ do i=1,ntyp
+!c write (iout,*) "tu dochodze",i
+ read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
+ if (i.eq.10) then
+ dsc_inv(i)=0.0D0
+ else
+ dsc_inv(i)=1.0D0/dsc(i)
+ endif
+ if (i.ne.10) then
+ do j=1,nlob(i)
+ do k=1,3
+ do l=1,3
+ blower(l,k,j)=0.0D0
+ enddo
+ enddo
+ enddo
+ bsc(1,i)=0.0D0
+ read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
+ censc(1,1,-i)=censc(1,1,i)
+ censc(2,1,-i)=censc(2,1,i)
+ censc(3,1,-i)=-censc(3,1,i)
+ do j=2,nlob(i)
+ read (irotam,*) bsc(j,i)
+ read (irotam,*) (censc(k,j,i),k=1,3),&
+ ((blower(k,l,j),l=1,k),k=1,3)
+ censc(1,j,-i)=censc(1,j,i)
+ censc(2,j,-i)=censc(2,j,i)
+ censc(3,j,-i)=-censc(3,j,i)
+! BSC is amplitude of Gaussian
+ enddo
+ do j=1,nlob(i)
+ do k=1,3
+ do l=1,k
+ akl=0.0D0
+ do m=1,3
+ akl=akl+blower(k,m,j)*blower(l,m,j)
+ enddo
+ gaussc(k,l,j,i)=akl
+ gaussc(l,k,j,i)=akl
+ if (((k.eq.3).and.(l.ne.3)) &
+ .or.((l.eq.3).and.(k.ne.3))) then
+ gaussc(k,l,j,-i)=-akl
+ gaussc(l,k,j,-i)=-akl
+ else
+ gaussc(k,l,j,-i)=akl
+ gaussc(l,k,j,-i)=akl
+ endif
+ enddo
+ enddo
+ enddo
+ endif
+ enddo
+ close (irotam)
+ if (lprint) then
+ write (iout,'(/a)') 'Parameters of side-chain local geometry'
+ do i=1,ntyp
+ nlobi=nlob(i)
+ if (nlobi.gt.0) then
+ write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),&
+ ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
+! write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
+! write (iout,'(a,f10.4,4(16x,f10.4))')
+! & 'Center ',(bsc(j,i),j=1,nlobi)
+! write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
+ write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
+ 'log h',(bsc(j,i),j=1,nlobi)
+ write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
+ 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
+! write (iout,'(a)')
+! do j=1,nlobi
+! ind=0
+! do k=1,3
+! do l=1,k
+! ind=ind+1
+! blower(k,l,j)=gaussc(ind,j,i)
+! enddo
+! enddo
+! enddo
+ do k=1,3
+ write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
+ ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
+ enddo
+ endif
+ enddo
+ endif
+#else
+!
+! Read scrot parameters for potentials determined from all-atom AM1 calculations
+! added by Urszula Kozlowska 07/11/2007
+!
+ allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
+
+ do i=1,ntyp
+ read (irotam,*)
+ if (i.eq.10) then
+ read (irotam,*)
+ else
+ do j=1,65
+ read(irotam,*) sc_parmin(j,i)
+ enddo
+ endif
+ enddo
+#endif
+ close(irotam)
+#ifdef CRYST_TOR
+!
+! Read torsional parameters in old format
+!
+ allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
+
+ read (itorp,*) ntortyp,nterm_old
+ write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
+ read (itorp,*) (itortyp(i),i=1,ntyp)
+
+!el from energy module--------
+ allocate(v1(nterm_old,ntortyp,ntortyp))
+ allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
+!el---------------------------
+
+ do i=1,ntortyp
+ do j=1,ntortyp
+ read (itorp,'(a)')
+ do k=1,nterm_old
+ read (itorp,*) kk,v1(k,j,i),v2(k,j,i)
+ enddo
+ enddo
+ enddo
+ close (itorp)
+ if (lprint) then
+ write (iout,'(/a/)') 'Torsional constants:'
+ do i=1,ntortyp
+ do j=1,ntortyp
+ write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
+ write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
+ enddo
+ enddo
+ endif
+
+
+#else
+!
+! Read torsional parameters
+!
+ allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
+
+ read (itorp,*) ntortyp
+ read (itorp,*) (itortyp(i),i=1,ntyp)
+ write (iout,*) 'ntortyp',ntortyp
+
+!el from energy module---------
+ allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
+ allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
+
+ allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
+ allocate(vlor2(maxlor,ntortyp,ntortyp))
+ allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
+ allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
+
+ allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
+ allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
+!el---------------------------
+ do iblock=1,2
+ do i=-ntortyp,ntortyp
+ do j=-ntortyp,ntortyp
+ nterm(i,j,iblock)=0
+ nlor(i,j,iblock)=0
+ enddo
+ enddo
+ enddo
+!el---------------------------
+
+ do iblock=1,2
+ do i=-ntyp,-1
+ itortyp(i)=-itortyp(-i)
+ enddo
+! write (iout,*) 'ntortyp',ntortyp
+ do i=0,ntortyp-1
+ do j=-ntortyp+1,ntortyp-1
+ read (itorp,*) nterm(i,j,iblock),&
+ nlor(i,j,iblock)
+ nterm(-i,-j,iblock)=nterm(i,j,iblock)
+ nlor(-i,-j,iblock)=nlor(i,j,iblock)
+ v0ij=0.0d0
+ si=-1.0d0
+ do k=1,nterm(i,j,iblock)
+ read (itorp,*) kk,v1(k,i,j,iblock),&
+ v2(k,i,j,iblock)
+ v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
+ v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
+ v0ij=v0ij+si*v1(k,i,j,iblock)
+ si=-si
+ enddo
+ do k=1,nlor(i,j,iblock)
+ read (itorp,*) kk,vlor1(k,i,j),&
+ vlor2(k,i,j),vlor3(k,i,j)
+ v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
+ enddo
+ v0(i,j,iblock)=v0ij
+ v0(-i,-j,iblock)=v0ij
+ enddo
+ enddo
+ enddo
+ close (itorp)
+ if (lprint) then
+ do iblock=1,2 !el
+ write (iout,'(/a/)') 'Torsional constants:'
+ do i=1,ntortyp
+ do j=1,ntortyp
+ write (iout,*) 'ityp',i,' jtyp',j
+ write (iout,*) 'Fourier constants'
+ do k=1,nterm(i,j,iblock)
+ write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
+ v2(k,i,j,iblock)
+ enddo
+ write (iout,*) 'Lorenz constants'
+ do k=1,nlor(i,j,iblock)
+ write (iout,'(3(1pe15.5))') &
+ vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
+ enddo
+ enddo
+ enddo
+ enddo
+ endif
+!
+! 6/23/01 Read parameters for double torsionals
+!
+!el from energy module------------
+ allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
+ allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
+!(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
+ allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
+ allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
+ !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
+ allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
+ allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
+ !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
+!---------------------------------
+
+ do iblock=1,2
+ do i=0,ntortyp-1
+ do j=-ntortyp+1,ntortyp-1
+ do k=-ntortyp+1,ntortyp-1
+ read (itordp,'(3a1)') t1,t2,t3
+! write (iout,*) "OK onelett",
+! & i,j,k,t1,t2,t3
+
+ if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
+ .or. t3.ne.toronelet(k)) then
+ write (iout,*) "Error in double torsional parameter file",&
+ i,j,k,t1,t2,t3
+#ifdef MPI
+ call MPI_Finalize(Ierror)
+#endif
+ stop "Error in double torsional parameter file"
+ endif
+ read (itordp,*) ntermd_1(i,j,k,iblock),&
+ ntermd_2(i,j,k,iblock)
+ ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
+ ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
+ read (itordp,*) (v1c(1,l,i,j,k,iblock),l=1,&
+ ntermd_1(i,j,k,iblock))
+ read (itordp,*) (v1s(1,l,i,j,k,iblock),l=1,&
+ ntermd_1(i,j,k,iblock))
+ read (itordp,*) (v1c(2,l,i,j,k,iblock),l=1,&
+ ntermd_1(i,j,k,iblock))
+ read (itordp,*) (v1s(2,l,i,j,k,iblock),l=1,&
+ ntermd_1(i,j,k,iblock))
+! Martix of D parameters for one dimesional foureir series
+ do l=1,ntermd_1(i,j,k,iblock)
+ v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
+ v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
+ v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
+ v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
+! write(iout,*) "whcodze" ,
+! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
+ enddo
+ read (itordp,*) ((v2c(l,m,i,j,k,iblock),&
+ v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
+ v2s(m,l,i,j,k,iblock),&
+ m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
+! Martix of D parameters for two dimesional fourier series
+ do l=1,ntermd_2(i,j,k,iblock)
+ do m=1,l-1
+ v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
+ v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
+ v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
+ v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
+ enddo!m
+ enddo!l
+ enddo!k
+ enddo!j
+ enddo!i
+ enddo!iblock
+ if (lprint) then
+ write (iout,*)
+ write (iout,*) 'Constants for double torsionals'
+ do iblock=1,2
+ do i=0,ntortyp-1
+ do j=-ntortyp+1,ntortyp-1
+ do k=-ntortyp+1,ntortyp-1
+ write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
+ ' nsingle',ntermd_1(i,j,k,iblock),&
+ ' ndouble',ntermd_2(i,j,k,iblock)
+ write (iout,*)
+ write (iout,*) 'Single angles:'
+ do l=1,ntermd_1(i,j,k,iblock)
+ write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
+ v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
+ v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
+ v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
+ enddo
+ write (iout,*)
+ write (iout,*) 'Pairs of angles:'
+ write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
+ do l=1,ntermd_2(i,j,k,iblock)
+ write (iout,'(i5,20f10.5)') &
+ l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
+ enddo
+ write (iout,*)
+ write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
+ do l=1,ntermd_2(i,j,k,iblock)
+ write (iout,'(i5,20f10.5)') &
+ l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
+ (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
+ enddo
+ write (iout,*)
+ enddo
+ enddo
+ enddo
+ enddo
+ endif
+#endif
+!elwrite(iout,*) "parmread kontrol sc-bb"
+! Read of Side-chain backbone correlation parameters
+! Modified 11 May 2012 by Adasko
+!CC
+!
+ read (isccor,*) nsccortyp
+
+ maxinter=3
+!c maxinter is maximum interaction sites
+!write(iout,*)"maxterm_sccor",maxterm_sccor
+!el from module energy-------------
+ allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
+ allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
+ allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
+ allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
+!-----------------------------------
+ allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
+!-----------------------------------
+ allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
+ allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
+ -nsccortyp:nsccortyp))
+ allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
+ -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
+ allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
+ -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
+!-----------------------------------
+ do i=-nsccortyp,nsccortyp
+ do j=-nsccortyp,nsccortyp
+ nterm_sccor(j,i)=0
+ enddo
+ enddo
+!-----------------------------------
+
+ read (isccor,*) (isccortyp(i),i=1,ntyp)
+ do i=-ntyp,-1
+ isccortyp(i)=-isccortyp(-i)
+ enddo
+ iscprol=isccortyp(20)
+! write (iout,*) 'ntortyp',ntortyp
+! maxinter=3
+!c maxinter is maximum interaction sites
+ do l=1,maxinter
+ do i=1,nsccortyp
+ do j=1,nsccortyp
+ read (isccor,*) &
+ nterm_sccor(i,j),nlor_sccor(i,j)
+ v0ijsccor=0.0d0
+ v0ijsccor1=0.0d0
+ v0ijsccor2=0.0d0
+ v0ijsccor3=0.0d0
+ si=-1.0d0
+ nterm_sccor(-i,j)=nterm_sccor(i,j)
+ nterm_sccor(-i,-j)=nterm_sccor(i,j)
+ nterm_sccor(i,-j)=nterm_sccor(i,j)
+ do k=1,nterm_sccor(i,j)
+ read (isccor,*) kk,v1sccor(k,l,i,j),&
+ v2sccor(k,l,i,j)
+ if (j.eq.iscprol) then
+ if (i.eq.isccortyp(10)) then
+ v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
+ v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
+ else
+ v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
+ +v2sccor(k,l,i,j)*dsqrt(0.75d0)
+ v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
+ +v1sccor(k,l,i,j)*dsqrt(0.75d0)
+ v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
+ v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
+ v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
+ v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
+ endif
+ else
+ if (i.eq.isccortyp(10)) then
+ v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
+ v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
+ else
+ if (j.eq.isccortyp(10)) then
+ v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
+ v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
+ else
+ v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
+ v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
+ v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
+ v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
+ v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
+ v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
+ endif
+ endif
+ endif
+ v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
+ v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
+ v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
+ v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
+ si=-si
+ enddo
+ do k=1,nlor_sccor(i,j)
+ read (isccor,*) kk,vlor1sccor(k,i,j),&
+ vlor2sccor(k,i,j),vlor3sccor(k,i,j)
+ v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
+ (1+vlor3sccor(k,i,j)**2)
+ enddo
+ v0sccor(l,i,j)=v0ijsccor
+ v0sccor(l,-i,j)=v0ijsccor1
+ v0sccor(l,i,-j)=v0ijsccor2
+ v0sccor(l,-i,-j)=v0ijsccor3
+ enddo
+ enddo
+ enddo
+ close (isccor)
+ if (lprint) then
+ write (iout,'(/a/)') 'Torsional constants of SCCORR:'
+ do i=1,nsccortyp
+ do j=1,nsccortyp
+ write (iout,*) 'ityp',i,' jtyp',j
+ write (iout,*) 'Fourier constants'
+ do k=1,nterm_sccor(i,j)
+ write (iout,'(2(1pe15.5))') &
+ (v1sccor(k,l,i,j),v2sccor(k,l,i,j),l=1,maxinter)
+ enddo
+ write (iout,*) 'Lorenz constants'
+ do k=1,nlor_sccor(i,j)
+ write (iout,'(3(1pe15.5))') &
+ vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
+ enddo
+ enddo
+ enddo
+ endif
+!
+! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
+! interaction energy of the Gly, Ala, and Pro prototypes.
+!
+ read (ifourier,*) nloctyp
+!el write(iout,*)"nloctyp",nloctyp
+!el from module energy-------
+ allocate(b1(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
+ allocate(b2(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
+ allocate(b1tilde(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
+ allocate(cc(2,2,-nloctyp-1:nloctyp+1))
+ allocate(dd(2,2,-nloctyp-1:nloctyp+1))
+ allocate(ee(2,2,-nloctyp-1:nloctyp+1))
+ allocate(ctilde(2,2,-nloctyp-1:nloctyp+1))
+ allocate(dtilde(2,2,-nloctyp-1:nloctyp+1)) !(2,2,-maxtor:maxtor)
+ do i=1,2
+ do ii=-nloctyp-1,nloctyp+1
+ b1(i,ii)=0.0d0
+ b2(i,ii)=0.0d0
+ b1tilde(i,ii)=0.0d0
+ do j=1,2
+ cc(j,i,ii)=0.0d0
+ dd(j,i,ii)=0.0d0
+ ee(j,i,ii)=0.0d0
+ ctilde(j,i,ii)=0.0d0
+ dtilde(j,i,ii)=0.0d0
+ enddo
+ enddo
+ enddo
+!--------------------------------
+ allocate(b(13,0:nloctyp))
+
+ do i=0,nloctyp-1
+ read (ifourier,*)
+ read (ifourier,*) (b(ii,i),ii=1,13)
+ if (lprint) then
+ write (iout,*) 'Type',i
+ write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
+ endif
+ B1(1,i) = b(3,i)
+ B1(2,i) = b(5,i)
+ B1(1,-i) = b(3,i)
+ B1(2,-i) = -b(5,i)
+! b1(1,i)=0.0d0
+! b1(2,i)=0.0d0
+ B1tilde(1,i) = b(3,i)
+ B1tilde(2,i) =-b(5,i)
+ B1tilde(1,-i) =-b(3,i)
+ B1tilde(2,-i) =b(5,i)
+! b1tilde(1,i)=0.0d0
+! b1tilde(2,i)=0.0d0
+ B2(1,i) = b(2,i)
+ B2(2,i) = b(4,i)
+ B2(1,-i) =b(2,i)
+ B2(2,-i) =-b(4,i)
+
+! b2(1,i)=0.0d0
+! b2(2,i)=0.0d0
+ CC(1,1,i)= b(7,i)
+ CC(2,2,i)=-b(7,i)
+ CC(2,1,i)= b(9,i)
+ CC(1,2,i)= b(9,i)
+ CC(1,1,-i)= b(7,i)
+ CC(2,2,-i)=-b(7,i)
+ CC(2,1,-i)=-b(9,i)
+ CC(1,2,-i)=-b(9,i)
+! CC(1,1,i)=0.0d0
+! CC(2,2,i)=0.0d0
+! CC(2,1,i)=0.0d0
+! CC(1,2,i)=0.0d0
+ Ctilde(1,1,i)=b(7,i)
+ Ctilde(1,2,i)=b(9,i)
+ Ctilde(2,1,i)=-b(9,i)
+ Ctilde(2,2,i)=b(7,i)
+ Ctilde(1,1,-i)=b(7,i)
+ Ctilde(1,2,-i)=-b(9,i)
+ Ctilde(2,1,-i)=b(9,i)
+ Ctilde(2,2,-i)=b(7,i)
+
+! Ctilde(1,1,i)=0.0d0
+! Ctilde(1,2,i)=0.0d0
+! Ctilde(2,1,i)=0.0d0
+! Ctilde(2,2,i)=0.0d0
+ DD(1,1,i)= b(6,i)
+ DD(2,2,i)=-b(6,i)
+ DD(2,1,i)= b(8,i)
+ DD(1,2,i)= b(8,i)
+ DD(1,1,-i)= b(6,i)
+ DD(2,2,-i)=-b(6,i)
+ DD(2,1,-i)=-b(8,i)
+ DD(1,2,-i)=-b(8,i)
+! DD(1,1,i)=0.0d0
+! DD(2,2,i)=0.0d0
+! DD(2,1,i)=0.0d0
+! DD(1,2,i)=0.0d0
+ Dtilde(1,1,i)=b(6,i)
+ Dtilde(1,2,i)=b(8,i)
+ Dtilde(2,1,i)=-b(8,i)
+ Dtilde(2,2,i)=b(6,i)
+ Dtilde(1,1,-i)=b(6,i)
+ Dtilde(1,2,-i)=-b(8,i)
+ Dtilde(2,1,-i)=b(8,i)
+ Dtilde(2,2,-i)=b(6,i)
+
+! Dtilde(1,1,i)=0.0d0
+! Dtilde(1,2,i)=0.0d0
+! Dtilde(2,1,i)=0.0d0
+! Dtilde(2,2,i)=0.0d0
+ EE(1,1,i)= b(10,i)+b(11,i)
+ EE(2,2,i)=-b(10,i)+b(11,i)
+ EE(2,1,i)= b(12,i)-b(13,i)
+ EE(1,2,i)= b(12,i)+b(13,i)
+ EE(1,1,-i)= b(10,i)+b(11,i)
+ EE(2,2,-i)=-b(10,i)+b(11,i)
+ EE(2,1,-i)=-b(12,i)+b(13,i)
+ EE(1,2,-i)=-b(12,i)-b(13,i)
+
+! ee(1,1,i)=1.0d0
+! ee(2,2,i)=1.0d0
+! ee(2,1,i)=0.0d0
+! ee(1,2,i)=0.0d0
+! ee(2,1,i)=ee(1,2,i)
+
+ enddo
+ if (lprint) then
+ do i=1,nloctyp
+ write (iout,*) 'Type',i
+ write (iout,*) 'B1'
+! write (iout,'(f10.5)') B1(:,i)
+ write(iout,*) B1(1,i),B1(2,i)
+ write (iout,*) 'B2'
+! write (iout,'(f10.5)') B2(:,i)
+ write(iout,*) B2(1,i),B2(2,i)
+ write (iout,*) 'CC'
+ do j=1,2
+ write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
+ enddo
+ write(iout,*) 'DD'
+ do j=1,2
+ write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
+ enddo
+ write(iout,*) 'EE'
+ do j=1,2
+ write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
+ enddo
+ enddo
+ endif
+!
+! Read electrostatic-interaction parameters
+!
+ if (lprint) then
+ write (iout,'(/a)') 'Electrostatic interaction constants:'
+ write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
+ 'IT','JT','APP','BPP','AEL6','AEL3'
+ endif
+ read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
+ read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
+ read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
+ read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
+ close (ielep)
+ do i=1,2
+ do j=1,2
+ rri=rpp(i,j)**6
+ app (i,j)=epp(i,j)*rri*rri
+ bpp (i,j)=-2.0D0*epp(i,j)*rri
+ ael6(i,j)=elpp6(i,j)*4.2D0**6
+ ael3(i,j)=elpp3(i,j)*4.2D0**3
+ if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
+ ael6(i,j),ael3(i,j)
+ enddo
+ enddo
+!
+! Read side-chain interaction parameters.
+!
+!el from module energy - COMMON.INTERACT-------
+ allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
+ allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
+ allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
+ allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
+ allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
+ do i=1,ntyp
+ do j=1,ntyp
+ augm(i,j)=0.0D0
+ enddo
+ chip(i)=0.0D0
+ alp(i)=0.0D0
+ sigma0(i)=0.0D0
+ sigii(i)=0.0D0
+ rr0(i)=0.0D0
+ enddo
+!--------------------------------
+
+ read (isidep,*) ipot,expon
+!el if (ipot.lt.1 .or. ipot.gt.5) then
+! write (iout,'(2a)') 'Error while reading SC interaction',&
+! 'potential file - unknown potential type.'
+! stop
+!wl endif
+ expon2=expon/2
+ write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
+ ', exponents are ',expon,2*expon
+! goto (10,20,30,30,40) ipot
+ select case(ipot)
+!----------------------- LJ potential ---------------------------------
+ case (1)
+! 10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
+ read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
+ if (lprint) then
+ write (iout,'(/a/)') 'Parameters of the LJ potential:'
+ write (iout,'(a/)') 'The epsilon array:'
+ call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
+ write (iout,'(/a)') 'One-body parameters:'
+ write (iout,'(a,4x,a)') 'residue','sigma'
+ write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
+ endif
+! goto 50
+!----------------------- LJK potential --------------------------------
+ case (2)
+! 20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
+ read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
+ (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
+ if (lprint) then
+ write (iout,'(/a/)') 'Parameters of the LJK potential:'
+ write (iout,'(a/)') 'The epsilon array:'
+ call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
+ write (iout,'(/a)') 'One-body parameters:'
+ write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
+ write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),&
+ i=1,ntyp)
+ endif
+! goto 50
+!---------------------- GB or BP potential -----------------------------
+ case (3:4)
+! 30 do i=1,ntyp
+ do i=1,ntyp
+ read (isidep,*)(eps(i,j),j=i,ntyp)
+ enddo
+ read (isidep,*)(sigma0(i),i=1,ntyp)
+ read (isidep,*)(sigii(i),i=1,ntyp)
+ read (isidep,*)(chip(i),i=1,ntyp)
+ read (isidep,*)(alp(i),i=1,ntyp)
+! For the GB potential convert sigma'**2 into chi'
+ if (ipot.eq.4) then
+ do i=1,ntyp
+ chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
+ enddo
+ endif
+ if (lprint) then
+ write (iout,'(/a/)') 'Parameters of the BP potential:'
+ write (iout,'(a/)') 'The epsilon array:'
+ call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
+ write (iout,'(/a)') 'One-body parameters:'
+ write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
+ ' chip ',' alph '
+ write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),&
+ chip(i),alp(i),i=1,ntyp)
+ endif
+! goto 50
+!--------------------- GBV potential -----------------------------------
+ case (5)
+! 40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
+ read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),&
+ (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
+ (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
+ if (lprint) then
+ write (iout,'(/a/)') 'Parameters of the GBV potential:'
+ write (iout,'(a/)') 'The epsilon array:'
+ call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
+ write (iout,'(/a)') 'One-body parameters:'
+ write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
+ 's||/s_|_^2',' chip ',' alph '
+ write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),&
+ sigii(i),chip(i),alp(i),i=1,ntyp)
+ endif
+ case default
+ write (iout,'(2a)') 'Error while reading SC interaction',&
+ 'potential file - unknown potential type.'
+ stop
+! 50 continue
+ end select
+! continue
+ close (isidep)
+!-----------------------------------------------------------------------
+! Calculate the "working" parameters of SC interactions.
+
+!el from module energy - COMMON.INTERACT-------
+ allocate(aa(ntyp1,ntyp1),bb(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
+ allocate(sigma(0:ntyp1,0:ntyp1),r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
+ do i=1,ntyp1
+ do j=1,ntyp1
+ aa(i,j)=0.0D0
+ bb(i,j)=0.0D0
+ chi(i,j)=0.0D0
+ sigma(i,j)=0.0D0
+ r0(i,j)=0.0D0
+ enddo
+ enddo
+!--------------------------------
+
+ do i=2,ntyp
+ do j=1,i-1
+ eps(i,j)=eps(j,i)
+ enddo
+ enddo
+ do i=1,ntyp
+ do j=i,ntyp
+ sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
+ sigma(j,i)=sigma(i,j)
+ rs0(i,j)=dwa16*sigma(i,j)
+ rs0(j,i)=rs0(i,j)
+ enddo
+ enddo
+ if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
+ 'Working parameters of the SC interactions:',&
+ ' a ',' b ',' augm ',' sigma ',' r0 ',&
+ ' chi1 ',' chi2 '
+ do i=1,ntyp
+ do j=i,ntyp
+ epsij=eps(i,j)
+ if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
+ rrij=sigma(i,j)
+ else
+ rrij=rr0(i)+rr0(j)
+ endif
+ r0(i,j)=rrij
+ r0(j,i)=rrij
+ rrij=rrij**expon
+ epsij=eps(i,j)
+ sigeps=dsign(1.0D0,epsij)
+ epsij=dabs(epsij)
+ aa(i,j)=epsij*rrij*rrij
+ bb(i,j)=-sigeps*epsij*rrij
+ aa(j,i)=aa(i,j)
+ bb(j,i)=bb(i,j)
+ if (ipot.gt.2) then
+ sigt1sq=sigma0(i)**2
+ sigt2sq=sigma0(j)**2
+ sigii1=sigii(i)
+ sigii2=sigii(j)
+ ratsig1=sigt2sq/sigt1sq
+ ratsig2=1.0D0/ratsig1
+ chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
+ if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
+ rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
+ else
+ rsum_max=sigma(i,j)
+ endif
+! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
+ sigmaii(i,j)=rsum_max
+ sigmaii(j,i)=rsum_max
+! else
+! sigmaii(i,j)=r0(i,j)
+! sigmaii(j,i)=r0(i,j)
+! endif
+!d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
+ if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
+ r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
+ augm(i,j)=epsij*r_augm**(2*expon)
+! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
+ augm(j,i)=augm(i,j)
+ else
+ augm(i,j)=0.0D0
+ augm(j,i)=0.0D0
+ endif
+ if (lprint) then
+ write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
+ restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),&
+ sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
+ endif
+ enddo
+ enddo
+
+ allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
+ do i=1,ntyp
+ do j=1,2
+ bad(i,j)=0.0D0
+ enddo
+ enddo
+#ifdef CLUSTER
+!
+! Define the SC-p interaction constants
+!
+ do i=1,20
+ do j=1,2
+ eps_scp(i,j)=-1.5d0
+ rscp(i,j)=4.0d0
+ enddo
+ enddo
+#endif
+
+!elwrite(iout,*) "parmread kontrol before oldscp"
+!
+! Define the SC-p interaction constants
+!
+#ifdef OLDSCP
+ do i=1,20
+! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
+! helix formation)
+! aad(i,1)=0.3D0*4.0D0**12
+! Following line for constants currently implemented
+! "Hard" SC-p repulsion (gives correct turn spacing in helices)
+ aad(i,1)=1.5D0*4.0D0**12
+! aad(i,1)=0.17D0*5.6D0**12
+ aad(i,2)=aad(i,1)
+! "Soft" SC-p repulsion
+ bad(i,1)=0.0D0
+! Following line for constants currently implemented
+! aad(i,1)=0.3D0*4.0D0**6
+! "Hard" SC-p repulsion
+ bad(i,1)=3.0D0*4.0D0**6
+! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
+ bad(i,2)=bad(i,1)
+! aad(i,1)=0.0D0
+! aad(i,2)=0.0D0
+! bad(i,1)=1228.8D0
+! bad(i,2)=1228.8D0
+ enddo
+#else
+!
+! 8/9/01 Read the SC-p interaction constants from file
+!
+ do i=1,ntyp
+ read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
+ enddo
+ do i=1,ntyp
+ aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
+ aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
+ bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
+ bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
+ enddo
+
+ if (lprint) then
+ write (iout,*) "Parameters of SC-p interactions:"
+ do i=1,20
+ write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
+ eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
+ enddo
+ endif
+#endif
+!
+! Define the constants of the disulfide bridge
+!
+ ebr=-5.50D0
+!
+! Old arbitrary potential - commented out.
+!
+! dbr= 4.20D0
+! fbr= 3.30D0
+!
+! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
+! energy surface of diethyl disulfide.
+! A. Liwo and U. Kozlowska, 11/24/03
+!
+ D0CM = 3.78d0
+ AKCM = 15.1d0
+ AKTH = 11.0d0
+ AKCT = 12.0d0
+ V1SS =-1.08d0
+ V2SS = 7.61d0
+ V3SS = 13.7d0
+
+ if (lprint) then
+ write (iout,'(/a)') "Disulfide bridge parameters:"
+ write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
+ write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
+ write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
+ write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
+ ' v3ss:',v3ss
+ endif
+ return
+ end subroutine parmread
+#ifndef CLUSTER
+!-----------------------------------------------------------------------------
+! mygetenv.F
+!-----------------------------------------------------------------------------
+ subroutine mygetenv(string,var)
+!
+! Version 1.0
+!
+! This subroutine passes the environmental variables to FORTRAN program.
+! If the flags -DMYGETENV and -DMPI are not for compilation, it calls the
+! standard FORTRAN GETENV subroutine. If both flags are set, the subroutine
+! reads the environmental variables from $HOME/.env
+!
+! Usage: As for the standard FORTRAN GETENV subroutine.
+!
+! Purpose: some versions/installations of MPI do not transfer the environmental
+! variables to slave processors, if these variables are set in the shell script
+! from which mpirun is called.
+!
+! A.Liwo, 7/29/01
+!
+#ifdef MPI
+ use MPI_data
+ include "mpif.h"
+#endif
+! implicit none
+ character*(*) :: string,var
+#if defined(MYGETENV) && defined(MPI)
+! include "DIMENSIONS.ZSCOPT"
+! include "mpif.h"
+! include "COMMON.MPI"
+!el character*360 ucase
+!el external ucase
+ character(len=360) :: string1(360),karta
+ character(len=240) :: home
+ integer i,n !,ilen
+!el external ilen
+ call getenv("HOME",home)
+ open(99,file=home(:ilen(home))//"/.env",status="OLD",err=112)
+ do while (.true.)
+ read (99,end=111,err=111,'(a)') karta
+ do i=1,80
+ string1(i)=" "
+ enddo
+ call split_string(karta,string1,80,n)
+ if (ucase(string1(1)(:ilen(string1(1)))).eq."SETENV" .and. &
+ string1(2)(:ilen(string1(2))).eq.string(:ilen(string)) ) then
+ var=string1(3)
+ print *,"Processor",me,": ",var(:ilen(var)),&
+ " assigned to ",string(:ilen(string))
+ close(99)
+ return
+ endif
+ enddo
+ 111 print *,"Environment variable ",string(:ilen(string))," not set."
+ close(99)
+ return
+ 112 print *,"Error opening environment file!"
+#else
+ call getenv(string,var)
+#endif
+ return
+ end subroutine mygetenv
+!-----------------------------------------------------------------------------
+! readrtns.F
+!-----------------------------------------------------------------------------
+ subroutine read_general_data(*)
+
+ use control_data, only:indpdb,symetr
+ use energy_data, only:distchainmax
+! implicit none
+! include "DIMENSIONS"
+! include "DIMENSIONS.ZSCOPT"
+! include "DIMENSIONS.FREE"
+! include "COMMON.TORSION"
+! include "COMMON.INTERACT"
+! include "COMMON.IOUNITS"
+! include "COMMON.TIME1"
+! include "COMMON.PROT"
+! include "COMMON.PROTFILES"
+! include "COMMON.CHAIN"
+! include "COMMON.NAMES"
+! include "COMMON.FFIELD"
+! include "COMMON.ENEPS"
+! include "COMMON.WEIGHTS"
+! include "COMMON.FREE"
+! include "COMMON.CONTROL"
+! include "COMMON.ENERGIES"
+ character(len=800) :: controlcard
+ integer :: i,j,k,ii,n_ene_found
+ integer :: ind,itype1,itype2,itypf,itypsc,itypp
+!el integer ilen
+!el external ilen
+!el character*16 ucase
+ character(len=16) :: key
+!el external ucase
+ call card_concat(controlcard,.true.)
+ call readi(controlcard,"N_ENE",n_eneW,max_eneW)
+ if (n_eneW.gt.max_eneW) then
+ write (iout,*) "Error: parameter out of range: N_ENE",n_eneW,&
+ max_eneW
+ return 1
+ endif
+ call readi(controlcard,"NPARMSET",nparmset,1)
+write(iout,*)"in read_gen data"
+ separate_parset = index(controlcard,"SEPARATE_PARSET").gt.0
+ call readi(controlcard,"IPARMPRINT",iparmprint,1)
+ write (iout,*) "PARMPRINT",iparmprint
+ if (nparmset.gt.max_parm) then
+ write (iout,*) "Error: parameter out of range: NPARMSET",&
+ nparmset, Max_Parm
+ return 1
+ endif
+write(iout,*)"in read_gen data"
+ call readi(controlcard,"MAXIT",maxit,5000)
+ call reada(controlcard,"FIMIN",fimin,1.0d-3)
+ call readi(controlcard,"ENSEMBLES",ensembles,0)
+ hamil_rep=index(controlcard,"HAMIL_REP").gt.0
+ write (iout,*) "Number of energy parameter sets",nparmset
+ allocate(isampl(nparmset))
+ call multreadi(controlcard,"ISAMPL",isampl,nparmset,1)
+ write (iout,*) "MaxSlice",MaxSlice
+ call readi(controlcard,"NSLICE",nslice,1)
+write(iout,*)"in read_gen data"
+ call flush(iout)
+ if (nslice.gt.MaxSlice) then
+ write (iout,*) "Error: parameter out of range: NSLICE",nslice,&
+ MaxSlice
+ return 1
+ endif
+ write (iout,*) "Frequency of storing conformations",&
+ (isampl(i),i=1,nparmset)
+ write (iout,*) "Maxit",maxit," Fimin",fimin
+ call readi(controlcard,"NQ",nQ,1)
+ if (nQ.gt.MaxQ) then
+ write (iout,*) "Error: parameter out of range: NQ",nq,&
+ maxq
+ return 1
+ endif
+ indpdb=0
+ if (index(controlcard,"CLASSIFY").gt.0) indpdb=1
+ call reada(controlcard,"DELTA",delta,1.0d-2)
+ call readi(controlcard,"EINICHECK",einicheck,2)
+ call reada(controlcard,"DELTRMS",deltrms,5.0d-2)
+ call reada(controlcard,"DELTRGY",deltrgy,5.0d-2)
+ call readi(controlcard,"RESCALE",rescale_modeW,1)
+ check_conf=index(controlcard,"NO_CHECK_CONF").eq.0
+ call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
+ call readi(controlcard,'SYM',symetr,1)
+ write (iout,*) "DISTCHAINMAX",distchainmax
+ write (iout,*) "delta",delta
+ write (iout,*) "einicheck",einicheck
+ write (iout,*) "rescale_mode",rescale_modeW
+ call flush(iout)
+ bxfile=index(controlcard,"BXFILE").gt.0
+ cxfile=index(controlcard,"CXFILE").gt.0
+ if (nslice .eq. 1 .and. .not.bxfile .and. .not.cxfile) &
+ bxfile=.true.
+ histfile=index(controlcard,"HISTFILE").gt.0
+ histout=index(controlcard,"HISTOUT").gt.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)
+ return
+ end subroutine read_general_data
+!------------------------------------------------------------------------------
+ subroutine read_efree(*)
+!
+! Read molecular data
+!
+! implicit none
+! include 'DIMENSIONS'
+! include 'DIMENSIONS.ZSCOPT'
+! include 'DIMENSIONS.COMPAR'
+! include 'DIMENSIONS.FREE'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.TIME1'
+! include 'COMMON.SBRIDGE'
+! include 'COMMON.CONTROL'
+! include 'COMMON.CHAIN'
+! include 'COMMON.HEADER'
+! include 'COMMON.GEO'
+! include 'COMMON.FREE'
+ character(len=320) :: controlcard !,ucase
+ integer :: iparm,ib,i,j,npars
+!el integer ilen
+!el external ilen
+
+ if (hamil_rep) then
+ npars=1
+ else
+ npars=nParmSet
+ endif
+
+! call alloc_wham_arrays
+! allocate(nT_h(nParmSet))
+! allocate(replica(nParmSet))
+! allocate(umbrella(nParmSet))
+! allocate(read_iset(nParmSet))
+! allocate(nT_h(nParmSet))
+
+ do iparm=1,npars
+
+ call card_concat(controlcard,.true.)
+ call readi(controlcard,'NT',nT_h(iparm),1)
+ write (iout,*) "iparm",iparm," nt",nT_h(iparm)
+ call flush(iout)
+ if (nT_h(iparm).gt.MaxT_h) then
+ write (iout,*) "Error: parameter out of range: NT",nT_h(iparm),&
+ MaxT_h
+ return 1
+ endif
+ replica(iparm)=index(controlcard,"REPLICA").gt.0
+ umbrella(iparm)=index(controlcard,"UMBRELLA").gt.0
+ read_iset(iparm)=index(controlcard,"READ_ISET").gt.0
+ write (iout,*) "nQ",nQ," nT",nT_h(iparm)," replica ",&
+ replica(iparm)," umbrella ",umbrella(iparm),&
+ " read_iset",read_iset(iparm)
+ call flush(iout)
+ do ib=1,nT_h(iparm)
+ call card_concat(controlcard,.true.)
+ call readi(controlcard,'NR',nR(ib,iparm),1)
+ if (umbrella(iparm)) then
+ nRR(ib,iparm)=1
+ else
+ nRR(ib,iparm)=nR(ib,iparm)
+ endif
+ if (nR(ib,iparm).gt.MaxR) then
+ write (iout,*) "Error: parameter out of range: NR",&
+ nR(ib,iparm),MaxR
+ return 1
+ endif
+ call reada(controlcard,'TEMP',beta_h(ib,iparm),298.0d0)
+ beta_h(ib,iparm)=1.0d0/(beta_h(ib,iparm)*1.987D-3)
+ call multreada(controlcard,'FI',f(1,ib,iparm),nR(ib,iparm),&
+ 0.0d0)
+ do i=1,nR(ib,iparm)
+ call card_concat(controlcard,.true.)
+ call multreada(controlcard,'KH',KH(1,i,ib,iparm),nQ,&
+ 100.0d0)
+ call multreada(controlcard,'Q0',Q0(1,i,ib,iparm),nQ,&
+ 0.0d0)
+ enddo
+ enddo
+ do ib=1,nT_h(iparm)
+ write (iout,*) "ib",ib," beta_h",&
+ 1.0d0/(0.001987*beta_h(ib,iparm))
+ write (iout,*) "nR",nR(ib,iparm)
+ write (iout,*) "fi",(f(i,ib,iparm),i=1,nR(ib,iparm))
+ do i=1,nR(ib,iparm)
+ write (iout,*) "i",i," Kh",(Kh(j,i,ib,iparm),j=1,nQ),&
+ "q0",(q0(j,i,ib,iparm),j=1,nQ)
+ enddo
+ call flush(iout)
+ enddo
+
+ enddo
+
+ if (hamil_rep) then
+
+ do iparm=2,nParmSet
+ nT_h(iparm)=nT_h(1)
+ do ib=1,nT_h(iparm)
+ nR(ib,iparm)=nR(ib,1)
+ if (umbrella(iparm)) then
+ nRR(ib,iparm)=1
+ else
+ nRR(ib,iparm)=nR(ib,1)
+ endif
+ beta_h(ib,iparm)=beta_h(ib,1)
+ do i=1,nR(ib,iparm)
+ f(i,ib,iparm)=f(i,ib,1)
+ do j=1,nQ
+ KH(j,i,ib,iparm)=KH(j,i,ib,1)
+ Q0(j,i,ib,iparm)=Q0(j,i,ib,1)
+ enddo
+ enddo
+ replica(iparm)=replica(1)
+ umbrella(iparm)=umbrella(1)
+ read_iset(iparm)=read_iset(1)
+ enddo
+ enddo
+
+ endif
+
+ return
+ end subroutine read_efree
+!-----------------------------------------------------------------------------
+ subroutine read_protein_data(*)
+! implicit none
+! include "DIMENSIONS"
+! include "DIMENSIONS.ZSCOPT"
+! include "DIMENSIONS.FREE"
+#ifdef MPI
+ use MPI_data
+ include "mpif.h"
+ integer :: IERROR,ERRCODE!,STATUS(MPI_STATUS_SIZE)
+! include "COMMON.MPI"
+#endif
+! include "COMMON.CHAIN"
+! include "COMMON.IOUNITS"
+! include "COMMON.PROT"
+! include "COMMON.PROTFILES"
+! include "COMMON.NAMES"
+! include "COMMON.FREE"
+! include "COMMON.OBCINKA"
+ character(len=64) :: nazwa
+ character(len=16000) :: controlcard
+ integer :: i,ii,ib,iR,iparm,nthr,npars !,ilen,iroof
+!el external ilen,iroof
+ if (hamil_rep) then
+ npars=1
+ else
+ npars=nparmset
+ endif
+
+ do iparm=1,npars
+
+! Read names of files with conformation data.
+ if (replica(iparm)) then
+ nthr = 1
+ else
+ nthr = nT_h(iparm)
+ endif
+ do ib=1,nthr
+ do ii=1,nRR(ib,iparm)
+ write (iout,*) "Parameter set",iparm," temperature",ib,&
+ " window",ii
+ call flush(iout)
+ call card_concat(controlcard,.true.)
+ write (iout,*) controlcard(:ilen(controlcard))
+ call readi(controlcard,"NFILE_BIN",nfile_bin(ii,ib,iparm),0)
+ call readi(controlcard,"NFILE_ASC",nfile_asc(ii,ib,iparm),0)
+ call readi(controlcard,"NFILE_CX",nfile_cx(ii,ib,iparm),0)
+ call readi(controlcard,"REC_START",rec_start(ii,ib,iparm),1)
+ call readi(controlcard,"REC_END",rec_end(ii,ib,iparm),&
+ maxstr*isampl(iparm)+rec_start(ii,ib,iparm)-1)
+ call reada(controlcard,"TIME_START",&
+ time_start_collect(ii,ib,iparm),0.0d0)
+ call reada(controlcard,"TIME_END",time_end_collect(ii,ib,iparm),&
+ 1.0d10)
+ write (iout,*) "rec_start",rec_start(ii,ib,iparm),&
+ " rec_end",rec_end(ii,ib,iparm)
+ write (iout,*) "time_start",time_start_collect(ii,ib,iparm),&
+ " time_end",time_end_collect(ii,ib,iparm)
+ call flush(iout)
+ if (replica(iparm)) then
+ call readi(controlcard,"TOTRAJ",totraj(ii,iparm),1)
+ write (iout,*) "Number of trajectories",totraj(ii,iparm)
+ call flush(iout)
+ endif
+ if (nfile_bin(ii,ib,iparm).lt.2 &
+ .and. nfile_asc(ii,ib,iparm).eq.0 &
+ .and. nfile_cx(ii,ib,iparm).eq.0) then
+ write (iout,*) "Error - no action specified!"
+ return 1
+ endif
+ if (nfile_bin(ii,ib,iparm).gt.0) then
+ call card_concat(controlcard,.false.)
+ call split_string(controlcard,protfiles(1,1,ii,ib,iparm),&
+ maxfile_prot,nfile_bin(ii,ib,iparm))
+#ifdef DEBUG
+ write(iout,*)"nfile_bin",nfile_bin(ii,ib,iparm)
+ write(iout,*) (protfiles(i,1,ii,ib,iparm),&
+ i=1,nfile_bin(ii,ib,iparm))
+#endif
+ endif
+ if (nfile_asc(ii,ib,iparm).gt.0) then
+ call card_concat(controlcard,.false.)
+ call split_string(controlcard,protfiles(1,2,ii,ib,iparm),&
+ maxfile_prot,nfile_asc(ii,ib,iparm))
+#ifdef DEBUG
+ write(iout,*) "nfile_asc(ii,ib,iparm)",nfile_asc(ii,ib,iparm)
+ write(iout,*) (protfiles(i,2,ii,ib,iparm),&
+ i=1,nfile_asc(ii,ib,iparm))
+#endif
+ else if (nfile_cx(ii,ib,iparm).gt.0) then
+ call card_concat(controlcard,.false.)
+ call split_string(controlcard,protfiles(1,2,ii,ib,iparm),&
+ maxfile_prot,nfile_cx(ii,ib,iparm))
+#ifdef DEBUG
+ write(iout,*) "nfile_cx(ii,ib,iparm)",nfile_cx(ii,ib,iparm)
+ write(iout,*) (protfiles(i,2,ii,ib,iparm),&
+ i=1,nfile_cx(ii,ib,iparm))
+#endif
+ endif
+ call flush(iout)
+ enddo
+ enddo
+
+ enddo
+ return
+ end subroutine read_protein_data
+!-------------------------------------------------------------------------------
+ subroutine readsss(rekord,lancuch,wartosc,default)
+! implicit none
+ character*(*) :: rekord,lancuch,wartosc,default
+ character(len=80) :: aux
+ integer :: lenlan,lenrec,iread,ireade
+!el external ilen
+!el logical iblnk
+!el external iblnk
+ lenlan=ilen(lancuch)
+ lenrec=ilen(rekord)
+ iread=index(rekord,lancuch(:lenlan)//"=")
+! print *,"rekord",rekord," lancuch",lancuch
+! print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
+ if (iread.eq.0) then
+ wartosc=default
+ return
+ endif
+ iread=iread+lenlan+1
+! print *,"iread",iread
+! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
+ do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
+ iread=iread+1
+! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
+ enddo
+! print *,"iread",iread
+ if (iread.gt.lenrec) then
+ wartosc=default
+ return
+ endif
+ ireade=iread+1
+! print *,"ireade",ireade
+ do while (ireade.lt.lenrec .and. &
+ .not.iblnk(rekord(ireade:ireade)))
+ ireade=ireade+1
+ enddo
+ wartosc=rekord(iread:ireade)
+ return
+ end subroutine readsss
+!----------------------------------------------------------------------------
+ subroutine multreads(rekord,lancuch,tablica,dim,default)
+! implicit none
+ integer :: dim,i
+ character*(*) rekord,lancuch,tablica(dim),default
+ character(len=80) :: aux
+ integer :: lenlan,lenrec,iread,ireade
+!el external ilen
+!el logical iblnk
+!el external iblnk
+ do i=1,dim
+ tablica(i)=default
+ enddo
+ lenlan=ilen(lancuch)
+ lenrec=ilen(rekord)
+ iread=index(rekord,lancuch(:lenlan)//"=")
+! print *,"rekord",rekord," lancuch",lancuch
+! print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
+ if (iread.eq.0) return
+ iread=iread+lenlan+1
+ do i=1,dim
+! print *,"iread",iread
+! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
+ do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
+ iread=iread+1
+! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
+ enddo
+! print *,"iread",iread
+ if (iread.gt.lenrec) return
+ ireade=iread+1
+! print *,"ireade",ireade
+ do while (ireade.lt.lenrec .and. &
+ .not.iblnk(rekord(ireade:ireade)))
+ ireade=ireade+1
+ enddo
+ tablica(i)=rekord(iread:ireade)
+ iread=ireade+1
+ enddo
+ end subroutine multreads
+!----------------------------------------------------------------------------
+ subroutine split_string(rekord,tablica,dim,nsub)
+! implicit none
+ integer :: dim,nsub,i,ii,ll,kk
+ character*(*) tablica(dim)
+ character*(*) rekord
+!el integer ilen
+!el external ilen
+ do i=1,dim
+ tablica(i)=" "
+ enddo
+ ii=1
+ ll = ilen(rekord)
+ nsub=0
+ do i=1,dim
+! Find the start of term name
+ kk = 0
+ do while (ii.le.ll .and. rekord(ii:ii).eq." ")
+ ii = ii+1
+ enddo
+! Parse the name into TABLICA(i) until blank found
+ do while (ii.le.ll .and. rekord(ii:ii).ne." ")
+ kk = kk+1
+ tablica(i)(kk:kk)=rekord(ii:ii)
+ ii = ii+1
+ enddo
+ if (kk.gt.0) nsub=nsub+1
+ if (ii.gt.ll) return
+ enddo
+ return
+ end subroutine split_string
+!--------------------------------------------------------------------------------
+! readrtns_compar.F
+!--------------------------------------------------------------------------------
+ subroutine read_compar
+!
+! Read molecular data
+!
+ use conform_compar, only:alloc_compar_arrays
+ use control_data, only:pdbref
+ use geometry_data, only:deg2rad,rad2deg
+! implicit none
+! include 'DIMENSIONS'
+! include 'DIMENSIONS.ZSCOPT'
+! include 'DIMENSIONS.COMPAR'
+! include 'DIMENSIONS.FREE'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.TIME1'
+! include 'COMMON.SBRIDGE'
+! include 'COMMON.CONTROL'
+! include 'COMMON.COMPAR'
+! include 'COMMON.CHAIN'
+! include 'COMMON.HEADER'
+! include 'COMMON.GEO'
+! include 'COMMON.FREE'
+ character(len=320) :: controlcard !,ucase
+ character(len=64) :: wfile
+!el integer ilen
+!el external ilen
+ integer :: i,j,k
+write(iout,*)"jestesmy w read_compar"
+ call card_concat(controlcard,.true.)
+ pdbref=(index(controlcard,'PDBREF').gt.0)
+ call reada(controlcard,'CUTOFF_UP',rmscut_base_up,4.0d0)
+ call reada(controlcard,'CUTOFF_LOW',rmscut_base_low,3.0d0)
+ call reada(controlcard,'RMSUP_LIM',rmsup_lim,4.0d0)
+ call reada(controlcard,'RMSUPUP_LIM',rmsupup_lim,7.5d0)
+ verbose = index(controlcard,"VERBOSE").gt.0
+ lgrp=index(controlcard,"STATIN").gt.0
+ lgrp_out=index(controlcard,"STATOUT").gt.0
+ merge_helices=index(controlcard,"DONT_MERGE_HELICES").eq.0
+ binary = index(controlcard,"BINARY").gt.0
+ rmscut_base_up=rmscut_base_up/50
+ rmscut_base_low=rmscut_base_low/50
+ call reada(controlcard,"FRAC_SEC",frac_sec,0.66666666d0)
+ call readi(controlcard,'NLEVEL',nlevel,1)
+ if (nlevel.lt.0) then
+ allocate(nfrag(2))
+ call alloc_compar_arrays(maxfrag,1)
+ goto 121
+ else
+ allocate(nfrag(nlevel))
+ endif
+! Read the data pertaining to elementary fragments (level 1)
+ call readi(controlcard,'NFRAG',nfrag(1),0)
+ write(iout,*)"nfrag(1)",nfrag(1)
+ call alloc_compar_arrays(nfrag(1),nlevel)
+ do j=1,nfrag(1)
+ call card_concat(controlcard,.true.)
+ write (iout,*) controlcard(:ilen(controlcard))
+ call readi(controlcard,'NPIECE',npiece(j,1),0)
+ call readi(controlcard,'N_SHIFT1',n_shift(1,j,1),0)
+ call readi(controlcard,'N_SHIFT2',n_shift(2,j,1),0)
+ call reada(controlcard,'ANGCUT',ang_cut(j),50.0d0)
+ call reada(controlcard,'MAXANG',ang_cut1(j),360.0d0)
+ call reada(controlcard,'FRAC_MIN',frac_min(j),0.666666d0)
+ call reada(controlcard,'NC_FRAC',nc_fragm(j,1),0.5d0)
+ call readi(controlcard,'NC_REQ',nc_req_setf(j,1),0)
+ call readi(controlcard,'RMS',irms(j,1),0)
+ call readi(controlcard,'LOCAL',iloc(j),1)
+ call readi(controlcard,'ELCONT',ielecont(j,1),1)
+ if (ielecont(j,1).eq.0) then
+ call readi(controlcard,'SCCONT',isccont(j,1),1)
+ endif
+ ang_cut(j)=ang_cut(j)*deg2rad
+ ang_cut1(j)=ang_cut1(j)*deg2rad
+ do k=1,npiece(j,1)
+ call card_concat(controlcard,.true.)
+ call readi(controlcard,'IFRAG1',ifrag(1,k,j),0)
+ call readi(controlcard,'IFRAG2',ifrag(2,k,j),0)
+ enddo
+ write(iout,*)"j",j," npiece",npiece(j,1)," ifrag",&
+ (ifrag(1,k,j),ifrag(2,k,j),&
+ k=1,npiece(j,1))," ang_cut",ang_cut(j)*rad2deg,&
+ " ang_cut1",ang_cut1(j)*rad2deg
+ write(iout,*)"n_shift",n_shift(1,j,1),n_shift(2,j,1)
+ write(iout,*)"nc_frac",nc_fragm(j,1)," nc_req",nc_req_setf(j,1)
+ write(iout,*)"irms",irms(j,1)," ielecont",ielecont(j,1),&
+ " ilocal",iloc(j)," isccont",isccont(j,1)
+ enddo
+! Read data pertaning to higher levels
+ do i=2,nlevel
+ call card_concat(controlcard,.true.)
+ call readi(controlcard,'NFRAG',NFRAG(i),0)
+ write (iout,*) "i",i," nfrag",nfrag(i)
+ do j=1,nfrag(i)
+ call card_concat(controlcard,.true.)
+ if (i.eq.2) then
+ call readi(controlcard,'ELCONT',ielecont(j,i),0)
+ if (ielecont(j,i).eq.0) then
+ call readi(controlcard,'SCCONT',isccont(j,i),1)
+ endif
+ call readi(controlcard,'RMS',irms(j,i),0)
+ else
+ ielecont(j,i)=0
+ isccont(j,i)=0
+ irms(j,i)=1
+ endif
+ call readi(controlcard,'NPIECE',npiece(j,i),0)
+ call readi(controlcard,'N_SHIFT1',n_shift(1,j,i),0)
+ call readi(controlcard,'N_SHIFT2',n_shift(2,j,i),0)
+ call multreadi(controlcard,'IPIECE',ipiece(1,j,i),&
+ npiece(j,i),0)
+ call reada(controlcard,'NC_FRAC',nc_fragm(j,i),0.5d0)
+ call readi(controlcard,'NC_REQ',nc_req_setf(j,i),0)
+ write(iout,*) "j",j," npiece",npiece(j,i)," n_shift",&
+ n_shift(1,j,i),n_shift(2,j,i)," ielecont",ielecont(j,i),&
+ " isccont",isccont(j,i)," irms",irms(j,i)
+ write(iout,*) "ipiece",(ipiece(k,j,i),k=1,npiece(j,i))
+ write(iout,*)"n_shift",n_shift(1,j,i),n_shift(2,j,i)
+ write(iout,*)"nc_frac",nc_fragm(j,i),&
+ " nc_req",nc_req_setf(j,i)
+ enddo
+ enddo
+ if (binary) write (iout,*) "Classes written in binary format."
+ return
+ 121 continue
+ call reada(controlcard,'ANGCUT_HEL',angcut_hel,50.0d0)
+ call reada(controlcard,'MAXANG_HEL',angcut1_hel,60.0d0)
+ call reada(controlcard,'ANGCUT_BET',angcut_bet,90.0d0)
+ call reada(controlcard,'MAXANG_BET',angcut1_bet,360.0d0)
+ call reada(controlcard,'ANGCUT_STRAND',angcut_strand,90.0d0)
+ call reada(controlcard,'MAXANG_STRAND',angcut1_strand,60.0d0)
+ call reada(controlcard,'FRAC_MIN',frac_min_set,0.666666d0)
+ call reada(controlcard,'NC_FRAC_HEL',ncfrac_hel,0.5d0)
+ call readi(controlcard,'NC_REQ_HEL',ncreq_hel,0)
+ call reada(controlcard,'NC_FRAC_BET',ncfrac_bet,0.5d0)
+ call reada(controlcard,'NC_FRAC_PAIR',ncfrac_pair,0.3d0)
+ call readi(controlcard,'NC_REQ_BET',ncreq_bet,0)
+ call readi(controlcard,'NC_REQ_PAIR',ncreq_pair,0)
+ call readi(controlcard,'NSHIFT_HEL',nshift_hel,3)
+ call readi(controlcard,'NSHIFT_BET',nshift_bet,3)
+ call readi(controlcard,'NSHIFT_STRAND',nshift_strand,3)
+ call readi(controlcard,'NSHIFT_PAIR',nshift_pair,3)
+ call readi(controlcard,'RMS_SINGLE',irms_single,0)
+ call readi(controlcard,'CONT_SINGLE',icont_single,1)
+ call readi(controlcard,'LOCAL_SINGLE',iloc_single,1)
+ call readi(controlcard,'RMS_PAIR',irms_pair,0)
+ call readi(controlcard,'CONT_PAIR',icont_pair,1)
+ call readi(controlcard,'SPLIT_BET',isplit_bet,0)
+ angcut_hel=angcut_hel*deg2rad
+ angcut1_hel=angcut1_hel*deg2rad
+ angcut_bet=angcut_bet*deg2rad
+ angcut1_bet=angcut1_bet*deg2rad
+ angcut_strand=angcut_strand*deg2rad
+ angcut1_strand=angcut1_strand*deg2rad
+ write (iout,*) "Automatic detection of structural elements"
+ write (iout,*) 'NC_FRAC_HEL',ncfrac_hel,' NC_REQ_HEL',ncreq_hel,&
+ ' NC_FRAC_BET',ncfrac_bet,' NC_REQ_BET',ncreq_bet,&
+ ' RMS_SINGLE',irms_single,' CONT_SINGLE',icont_single,&
+ ' NC_FRAC_PAIR',ncfrac_pair,' NC_REQ_PAIR',ncreq_pair,&
+ ' RMS_PAIR',irms_pair,' CONT_PAIR',icont_pair,&
+ ' SPLIT_BET',isplit_bet
+ write (iout,*) 'NSHIFT_HEL',nshift_hel,' NSHIFT_BET',nshift_bet,&
+ ' NSHIFT_STRAND',nshift_strand,' NSHIFT_PAIR',nshift_pair
+ write (iout,*) 'ANGCUT_HEL',angcut_hel*rad2deg,&
+ ' MAXANG_HEL',angcut1_hel*rad2deg
+ write (iout,*) 'ANGCUT_BET',angcut_bet*rad2deg,&
+ ' MAXANG_BET',angcut1_bet*rad2deg
+ write (iout,*) 'ANGCUT_STRAND',angcut_strand*rad2deg,&
+ ' MAXANG_STRAND',angcut1_strand*rad2deg
+ write (iout,*) 'FRAC_MIN',frac_min_set
+ return
+ end subroutine read_compar
+!--------------------------------------------------------------------------------
+! read_ref_str.F
+!--------------------------------------------------------------------------------
+ subroutine read_ref_structure(*)
+!
+! Read the reference structure from the PDB file or from a PDB file or in the form of the dihedral
+! angles.
+!
+ use control_data, only:pdbref
+ use geometry_data, only:nres,cref,c,dc,nsup,dc_norm,nend_sup,&
+ nstart_sup,nstart_seq,nperm,nres0
+ use energy_data, only:nct,nnt,icont_ref,ncont_ref,itype
+ use compare, only:seq_comp !,contact,elecont
+ use geometry, only:chainbuild,dist
+ use io_config, only:readpdb
+!
+ use conform_compar, only:contact,elecont
+! implicit none
+! include 'DIMENSIONS'
+! include 'DIMENSIONS.ZSCOPT'
+! include 'DIMENSIONS.COMPAR'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.GEO'
+! include 'COMMON.VAR'
+! include 'COMMON.INTERACT'
+! include 'COMMON.LOCAL'
+! include 'COMMON.NAMES'
+! include 'COMMON.CHAIN'
+! include 'COMMON.FFIELD'
+! include 'COMMON.SBRIDGE'
+! include 'COMMON.HEADER'
+! include 'COMMON.CONTROL'
+! include 'COMMON.CONTACTS1'
+! include 'COMMON.PEPTCONT'
+! include 'COMMON.TIME1'
+! include 'COMMON.COMPAR'
+ character(len=4) :: sequence(nres)
+!el integer rescode
+!el real(kind=8) :: x(maxvar)
+ integer :: itype_pdb(nres)
+!el logical seq_comp
+ integer :: i,j,k,nres_pdb,iaux
+ real(kind=8) :: ddsc !el,dist
+ integer :: kkk !,ilen
+!el external ilen
+!
+ nres0=nres
+ write (iout,*) "pdbref",pdbref
+ if (pdbref) then
+ read(inp,'(a)') pdbfile
+ write (iout,'(2a,1h.)') 'PDB data will be read from file ',&
+ pdbfile(:ilen(pdbfile))
+ open(ipdbin,file=pdbfile,status='old',err=33)
+ goto 34
+ 33 write (iout,'(a)') 'Error opening PDB file.'
+ return 1
+ 34 continue
+ do i=1,nres
+ itype_pdb(i)=itype(i)
+ enddo
+write(iout,*)"jestesmy przed readpdb"
+ call readpdb
+ do i=1,nres
+ iaux=itype_pdb(i)
+ itype_pdb(i)=itype(i)
+ itype(i)=iaux
+ enddo
+ close (ipdbin)
+ do kkk=1,nperm
+ nres_pdb=nres
+ nres=nres0
+ nstart_seq=nnt
+ if (nsup.le.(nct-nnt+1)) then
+ do i=0,nct-nnt+1-nsup
+ if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),&
+ nsup)) then
+ do j=nnt+nsup-1,nnt,-1
+ do k=1,3
+ cref(k,nres+j+i,kkk)=cref(k,nres_pdb+j,kkk)
+ enddo
+ enddo
+ do j=nnt+nsup-1,nnt,-1
+ do k=1,3
+ cref(k,j+i,kkk)=cref(k,j,kkk)
+ enddo
+ phi_ref(j+i)=phi_ref(j)
+ theta_ref(j+i)=theta_ref(j)
+ alph_ref(j+i)=alph_ref(j)
+ omeg_ref(j+i)=omeg_ref(j)
+ enddo
+#ifdef DEBUG
+ do j=nnt,nct
+ write (iout,'(i5,3f10.5,5x,3f10.5)') &
+ j,(cref(k,j,kkk),k=1,3),(cref(k,j+nres,kkk),k=1,3)
+ enddo
+#endif
+ nstart_seq=nnt+i
+ nstart_sup=nnt+i
+ goto 111
+ endif
+ enddo
+ write (iout,'(a)') &
+ 'Error - sequences to be superposed do not match.'
+ return 1
+ else
+ do i=0,nsup-(nct-nnt+1)
+ if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),&
+ nct-nnt+1)) &
+ then
+ nstart_sup=nstart_sup+i
+ nsup=nct-nnt+1
+ goto 111
+ endif
+ enddo
+ write (iout,'(a)') &
+ 'Error - sequences to be superposed do not match.'
+ endif
+ enddo
+ 111 continue
+ write (iout,'(a,i5)') &
+ 'Experimental structure begins at residue',nstart_seq
+ else
+ call read_angles(inp,*38)
+ goto 39
+ 38 write (iout,'(a)') 'Error reading reference structure.'
+ return 1
+ 39 call chainbuild
+ kkk=1
+ nstart_sup=nnt
+ nstart_seq=nnt
+ nsup=nct-nnt+1
+ do i=1,2*nres
+ do j=1,3
+ cref(j,i,kkk)=c(j,i)
+ enddo
+ enddo
+ endif
+ nend_sup=nstart_sup+nsup-1
+ do i=1,2*nres
+ do j=1,3
+ c(j,i)=cref(j,i,kkk)
+ enddo
+ enddo
+ do i=1,nres
+ do j=1,3
+ dc(j,nres+i)=cref(j,nres+i,kkk)-cref(j,i,kkk)
+ enddo
+ if (itype(i).ne.10) then
+ ddsc = dist(i,nres+i)
+ do j=1,3
+ dc_norm(j,nres+i)=dc(j,nres+i)/ddsc
+ enddo
+ else
+ do j=1,3
+ dc_norm(j,nres+i)=0.0d0
+ enddo
+ endif
+! write (iout,*) "i",i," dc_norm",(dc_norm(k,nres+i),k=1,3),
+! " norm",dc_norm(1,nres+i)**2+dc_norm(2,nres+i)**2+
+! dc_norm(3,nres+i)**2
+ do j=1,3
+ dc(j,i)=c(j,i+1)-c(j,i)
+ enddo
+ ddsc = dist(i,i+1)
+ do j=1,3
+ dc_norm(j,i)=dc(j,i)/ddsc
+ enddo
+ enddo
+! print *,"Calling contact"
+ call contact(.true.,ncont_ref,icont_ref(1,1),&
+ nstart_sup,nend_sup)
+! print *,"Calling elecont"
+ call elecont(.true.,ncont_pept_ref,&
+ icont_pept_ref(1,1),&
+ nstart_sup,nend_sup)
+ write (iout,'(a,i3,a,i3,a,i3,a)') &
+ 'Number of residues to be superposed:',nsup,&
+ ' (from residue',nstart_sup,' to residue',&
+ nend_sup,').'
+ return
+ end subroutine read_ref_structure
+!--------------------------------------------------------------------------------
+! geomout.F
+!--------------------------------------------------------------------------------
+ subroutine pdboutW(ii,temp,efree,etot,entropy,rmsdev)
+
+ use geometry_data, only:nres,c
+ use energy_data, only:nss,nnt,nct,ihpb,jhpb,itype
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'DIMENSIONS.ZSCOPT'
+! include 'COMMON.CHAIN'
+! include 'COMMON.INTERACT'
+! include 'COMMON.NAMES'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.HEADER'
+! include 'COMMON.SBRIDGE'
+ character(len=50) :: tytul
+ character(len=1),dimension(10) :: chainid=reshape((/'A','B','C',&
+ 'D','E','F','G','H','I','J'/),shape(chainid))
+ integer,dimension(nres) :: ica !(maxres)
+ real(kind=8) :: temp,efree,etot,entropy,rmsdev
+ integer :: ii,i,j,iti,ires,iatom,ichain
+ write(ipdb,'("REMARK CONF",i8," TEMPERATURE",f7.1," RMS",0pf7.2)')&
+ ii,temp,rmsdev
+ write (ipdb,'("REMARK DIMENSIONLESS FREE ENERGY",1pe15.5)') &
+ efree
+ write (ipdb,'("REMARK ENERGY",1pe15.5," ENTROPY",1pe15.5)') &
+ etot,entropy
+ iatom=0
+ ichain=1
+ ires=0
+ do i=nnt,nct
+ iti=itype(i)
+ if (iti.eq.ntyp1) then
+ ichain=ichain+1
+ ires=0
+ write (ipdb,'(a)') 'TER'
+ else
+ ires=ires+1
+ iatom=iatom+1
+ ica(i)=iatom
+ write (ipdb,10) iatom,restyp(iti),chainid(ichain),&
+ ires,(c(j,i),j=1,3)
+ if (iti.ne.10) then
+ iatom=iatom+1
+ write (ipdb,20) iatom,restyp(iti),chainid(ichain),&
+ ires,(c(j,nres+i),j=1,3)
+ endif
+ endif
+ enddo
+ write (ipdb,'(a)') 'TER'
+ do i=nnt,nct-1
+ if (itype(i).eq.ntyp1) cycle
+ if (itype(i).eq.10 .and. itype(i+1).ne.ntyp1) then
+ write (ipdb,30) ica(i),ica(i+1)
+ else if (itype(i).ne.10 .and. itype(i+1).ne.ntyp1) then
+ write (ipdb,30) ica(i),ica(i+1),ica(i)+1
+ else if (itype(i).ne.10 .and. itype(i+1).eq.ntyp1) then
+ write (ipdb,30) ica(i),ica(i)+1
+ endif
+ enddo
+ if (itype(nct).ne.10) then
+ write (ipdb,30) ica(nct),ica(nct)+1
+ endif
+ do i=1,nss
+ write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
+ enddo
+ write (ipdb,'(a6)') 'ENDMDL'
+ 10 FORMAT ('ATOM',I7,' CA ',A3,1X,A1,I4,4X,3F8.3,f15.3)
+ 20 FORMAT ('ATOM',I7,' CB ',A3,1X,A1,I4,4X,3F8.3,f15.3)
+ 30 FORMAT ('CONECT',8I5)
+ return
+ end subroutine pdboutW
+#endif
+!------------------------------------------------------------------------------
+ end module io_wham
+!-----------------------------------------------------------------------------
+!-----------------------------------------------------------------------------
+
projects / unres4.git / blobdiff