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: maxterm,maxlor,maxterm_sccor,& !maxtor 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(-nthetyp-1:nthetyp+1,& -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2)) !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2) allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,& -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2)) !(maxtheterm,-maxthetyp1:maxthetyp1,& ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2) allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,& -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2)) allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,& -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2)) allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,& -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2)) allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,& -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2)) !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,& ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2) allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,& -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2)) allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,& -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,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" aa0thet(:,:,:,:)=0.0d0 aathet(:,:,:,:,:)=0.0d0 bbthet(:,:,:,:,:,:)=0.0d0 ccthet(:,:,:,:,:,:)=0.0d0 ddthet(:,:,:,:,:,:)=0.0d0 eethet(:,:,:,:,:,:)=0.0d0 ffthet(:,:,:,:,:,:,:)=0.0d0 ggthet(:,:,:,:,:,:,:)=0.0d0 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) !elwrite(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 !elwrite(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) !elwrite(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 !elwrite(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 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 !----------------------------------------------------------------------------- !-----------------------------------------------------------------------------