+++ /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)
-!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
-!-----------------------------------------------------------------------------
-!-----------------------------------------------------------------------------
-
projects / unres4.git / blobdiff