--- /dev/null
+ module io_config
+
+ use names
+ use io_units
+ use io_base
+ use geometry_data
+ use geometry
+ use control_data, only:maxterm_sccor
+ implicit none
+!-----------------------------------------------------------------------------
+! Max. number of residue types and parameters in expressions for
+! virtual-bond angle bending potentials
+! integer,parameter :: maxthetyp=3
+! integer,parameter :: maxthetyp1=maxthetyp+1
+! ,maxtheterm=20,
+! & maxtheterm2=6,maxtheterm3=4,maxsingle=6,maxdouble=4,
+! & mmaxtheterm=maxtheterm)
+!-----------------------------------------------------------------------------
+! Max. number of types of dihedral angles & multiplicity of torsional barriers
+! and the number of terms in double torsionals
+! integer,parameter :: maxlor=3,maxtermd_1=8,maxtermd_2=8
+! parameter (maxtor=4,maxterm=10)
+!-----------------------------------------------------------------------------
+! Max number of torsional terms in SCCOR
+!el integer,parameter :: maxterm_sccor=6
+!-----------------------------------------------------------------------------
+ character(len=1),dimension(:),allocatable :: secstruc !(maxres)
+!-----------------------------------------------------------------------------
+!
+!
+!-----------------------------------------------------------------------------
+ contains
+#if !defined(WHAM_RUN) && !defined(CLUSTER)
+!-----------------------------------------------------------------------------
+! bank.F io_csa
+!-----------------------------------------------------------------------------
+ subroutine write_rbank(jlee,adif,nft)
+
+ use csa_data
+ use geometry_data, only: nres,rad2deg
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+! include 'COMMON.CHAIN'
+! include 'COMMON.GEO'
+!el local variables
+ integer :: nft,i,k,j,l,jlee
+ real(kind=8) :: adif
+
+ open(icsa_rbank,file=csa_rbank,status="unknown")
+ write (icsa_rbank,900) jlee,nbank,nstep,nft,icycle,adif
+ do k=1,nbank
+ write (icsa_rbank,952) k,rene(k),rrmsn(k),rpncn(k)
+ do j=1,numch
+ do l=2,nres-1
+ write (icsa_rbank,850) (rad2deg*rvar(i,l,j,k),i=1,4)
+ enddo
+ enddo
+ enddo
+ close(icsa_rbank)
+
+ 850 format (10f8.3)
+ 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
+ i8,i10,i2,f15.5)
+ 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
+ ' %NC ',0pf5.2)
+
+ return
+ end subroutine write_rbank
+!-----------------------------------------------------------------------------
+ subroutine read_rbank(jlee,adif)
+
+ use csa_data
+ use geometry_data, only: nres,deg2rad
+ use MPI_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+ include 'mpif.h'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+! include 'COMMON.CHAIN'
+! include 'COMMON.GEO'
+! include 'COMMON.SETUP'
+ character(len=80) :: karta
+!el local variables
+ integer :: nbankr,nstepr,nftr,icycler,kk,k,j,l,i,&
+ ierror,ierrcode,jlee,jleer
+ real(kind=8) :: adif
+
+ open(icsa_rbank,file=csa_rbank,status="old")
+ read (icsa_rbank,901) jleer,nbankr,nstepr,nftr,icycler,adif
+ print *,jleer,nbankr,nstepr,nftr,icycler,adif
+! print *, 'adif from read_rbank ',adif
+#ifdef MPI
+ if(nbankr.ne.nbank) then
+ write (iout,*) 'ERROR in READ_BANK: NBANKR',nbankr,&
+ ' NBANK',nbank
+ call mpi_abort(mpi_comm_world,ierror,ierrcode)
+ endif
+ if(jleer.ne.jlee) then
+ write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,&
+ ' JLEE',jlee
+ call mpi_abort(mpi_comm_world,ierror,ierrcode)
+ endif
+#endif
+
+ kk=0
+ do k=1,nbankr
+ read (icsa_rbank,'(a80)') karta
+ write(iout,*) "READ_RBANK: kk=",kk
+ write(iout,*) karta
+! if (index(karta,"*").gt.0) then
+! write (iout,*) "***** Stars in bankr ***** k=",k,
+! & " skipped"
+! do j=1,numch
+! do l=2,nres-1
+! read (30,850) (rdummy,i=1,4)
+! enddo
+! enddo
+! else
+ kk=kk+1
+ call reada(karta,"total E",rene(kk),1.0d20)
+ call reada(karta,"rmsd from N",rrmsn(kk),0.0d0)
+ call reada(karta,"%NC",rpncn(kk),0.0d0)
+ write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
+ "%NC",bpncn(kk),ibank(kk)
+! read (icsa_rbank,953) kdummy,rene(kk),rrmsn(kk),rpncn(kk)
+ do j=1,numch
+ do l=2,nres-1
+ read (icsa_rbank,850) (rvar(i,l,j,kk),i=1,4)
+! write (iout,850) (rvar(i,l,j,kk),i=1,4)
+ do i=1,4
+ rvar(i,l,j,kk)=deg2rad*rvar(i,l,j,kk)
+ enddo
+ enddo
+ enddo
+! endif
+ enddo
+!d write (*,*) "read_rbank ******************* kk",kk,
+!d & "nbankr",nbankr
+ if (kk.lt.nbankr) nbankr=kk
+!d do kk=1,nbankr
+!d print *,"kk=",kk
+!d do j=1,numch
+!d do l=2,nres-1
+!d write (*,850) (rvar(i,l,j,kk),i=1,4)
+!d enddo
+!d enddo
+!d enddo
+ close(icsa_rbank)
+
+ 850 format (10f8.3)
+ 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
+ 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2)
+
+ return
+ end subroutine read_rbank
+!-----------------------------------------------------------------------------
+ subroutine write_bank(jlee,nft)
+
+ use csa_data
+ use control_data, only: vdisulf
+ use geometry_data, only: nres,rad2deg
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+! include 'COMMON.CHAIN'
+! include 'COMMON.GEO'
+! include 'COMMON.SBRIDGE'
+! include 'COMMON.CONTROL'
+ character(len=7) :: chtmp
+ character(len=40) :: chfrm
+!el external ilen
+!el local variables
+ integer :: nft,k,l,i,j,jlee
+
+ open(icsa_bank,file=csa_bank,status="unknown")
+ write (icsa_bank,900) jlee,nbank,nstep,nft,icycle,cutdif
+ write (icsa_bank,902) nglob_csa, eglob_csa
+ open (igeom,file=intname,status='UNKNOWN')
+ do k=1,nbank
+ write (icsa_bank,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
+ if (vdisulf) write (icsa_bank,'(101i4)') &
+ bvar_nss(k),((bvar_ss(j,i,k),j=1,2),i=1,bvar_nss(k))
+ do j=1,numch
+ do l=2,nres-1
+ write (icsa_bank,850) (rad2deg*bvar(i,l,j,k),i=1,4)
+ enddo
+ enddo
+ if (bvar_nss(k).le.9) then
+ write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
+ bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
+ else
+ write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
+ bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
+ write (igeom,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
+ bvar_ss(2,i,k),i=10,bvar_nss(k))
+ endif
+ write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
+ write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
+ write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
+ write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
+ enddo
+ close(icsa_bank)
+ close(igeom)
+
+ if (nstep/200.gt.ilastnstep) then
+
+ ilastnstep=(ilastnstep+1)*1.5
+ write(chfrm,'(a2,i1,a1)') '(i',int(dlog10(dble(nstep))+1),')'
+ write(chtmp,chfrm) nstep
+ open(icsa_int,file=prefix(:ilen(prefix)) &
+ //'_'//chtmp(:ilen(chtmp))//'.int',status='UNKNOWN')
+ do k=1,nbank
+ if (bvar_nss(k).le.9) then
+ write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
+ bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
+ else
+ write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
+ bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
+ write (icsa_int,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
+ bvar_ss(2,i,k),i=10,bvar_nss(k))
+ endif
+ write (icsa_int,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
+ write (icsa_int,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
+ write (icsa_int,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
+ write (icsa_int,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
+ enddo
+ close(icsa_int)
+ endif
+
+
+ 200 format (8f10.4)
+ 850 format (10f8.3)
+ 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
+ i8,i10,i2,f15.5)
+ 902 format (1x,'nglob_csa =',i4,' eglob_csa =',1pe14.5)
+ 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
+ ' %NC ',0pf5.2,i5)
+
+ return
+ end subroutine write_bank
+!-----------------------------------------------------------------------------
+ subroutine write_bank_reminimized(jlee,nft)
+
+ use csa_data
+ use geometry_data, only: nres,rad2deg
+ use energy_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+! include 'COMMON.CHAIN'
+! include 'COMMON.GEO'
+! include 'COMMON.SBRIDGE'
+!el local variables
+ integer :: nft,i,l,j,k,jlee
+
+ open(icsa_bank_reminimized,file=csa_bank_reminimized,&
+ status="unknown")
+ write (icsa_bank_reminimized,900) &
+ jlee,nbank,nstep,nft,icycle,cutdif
+ open (igeom,file=intname,status='UNKNOWN')
+ do k=1,nbank
+ write (icsa_bank_reminimized,952) k,bene(k),brmsn(k),&
+ bpncn(k),ibank(k)
+ do j=1,numch
+ do l=2,nres-1
+ write (icsa_bank_reminimized,850) (rad2deg*bvar(i,l,j,k),i=1,4)
+ enddo
+ enddo
+ if (nss.le.9) then
+ write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
+ nss,(ihpb(i),jhpb(i),i=1,nss)
+ else
+ write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
+ nss,(ihpb(i),jhpb(i),i=1,9)
+ write (igeom,'(3X,11(1X,2I3))') (ihpb(i),jhpb(i),i=10,nss)
+ endif
+ write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
+ write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
+ write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
+ write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
+ enddo
+ close(icsa_bank_reminimized)
+ close(igeom)
+
+ 200 format (8f10.4)
+ 850 format (10f8.3)
+ 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
+ i8,i10,i2,f15.5)
+ 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
+ ' %NC ',0pf5.2,i5)
+
+ return
+ end subroutine write_bank_reminimized
+!-----------------------------------------------------------------------------
+ subroutine read_bank(jlee,nft,cutdifr)
+
+ use csa_data
+ use control_data, only: vdisulf
+ use geometry_data, only: nres,deg2rad
+ use energy_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+! include 'COMMON.CHAIN'
+! include 'COMMON.GEO'
+! include 'COMMON.CONTROL'
+! include 'COMMON.SBRIDGE'
+ character(len=80) :: karta
+! integer ilen
+!el external ilen
+!el local variables
+ integer :: nft,kk,k,l,i,j,jlee
+ real(kind=8) :: cutdifr
+
+ open(icsa_bank,file=csa_bank,status="old")
+ read (icsa_bank,901) jlee,nbank,nstep,nft,icycle,cutdifr
+ read (icsa_bank,902) nglob_csa, eglob_csa
+! if(jleer.ne.jlee) then
+! write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,
+! & ' JLEE',jlee
+! call mpi_abort(mpi_comm_world,ierror,ierrcode)
+! endif
+
+ kk=0
+ do k=1,nbank
+ read (icsa_bank,'(a80)') karta
+ write(iout,*) "READ_BANK: kk=",kk
+ write(iout,*) karta
+! if (index(karta,"*").gt.0) then
+! write (iout,*) "***** Stars in bank ***** k=",k,
+! & " skipped"
+! do j=1,numch
+! do l=2,nres-1
+! read (33,850) (rdummy,i=1,4)
+! enddo
+! enddo
+! else
+ kk=kk+1
+ call reada(karta,"total E",bene(kk),1.0d20)
+ call reada(karta,"rmsd from N",brmsn(kk),0.0d0)
+ call reada(karta,"%NC",bpncn(kk),0.0d0)
+ read (karta(ilen(karta)-1:),*,end=111,err=111) ibank(kk)
+ goto 112
+ 111 ibank(kk)=0
+ 112 continue
+ write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
+ "%NC",bpncn(kk),ibank(kk)
+! read (icsa_bank,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
+ if (vdisulf) then
+ read (icsa_bank,'(101i4)') &
+ bvar_nss(kk),((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
+ bvar_ns(kk)=ns-2*bvar_nss(kk)
+ write(iout,*) 'read SSBOND',bvar_nss(kk),&
+ ((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
+!d write(iout,*) 'read CYS #free ', bvar_ns(kk)
+ l=0
+ do i=1,ns
+ j=1
+ do while( iss(i).ne.bvar_ss(1,j,kk)-nres .and. &
+ iss(i).ne.bvar_ss(2,j,kk)-nres .and. &
+ j.le.bvar_nss(kk))
+ j=j+1
+ enddo
+ if (j.gt.bvar_nss(kk)) then
+ l=l+1
+ bvar_s(l,kk)=iss(i)
+ endif
+ enddo
+!d write(iout,*)'read CYS free',(bvar_s(l,kk),l=1,bvar_ns(kk))
+ endif
+ do j=1,numch
+ do l=2,nres-1
+ read (icsa_bank,850) (bvar(i,l,j,kk),i=1,4)
+! write (iout,850) (bvar(i,l,j,kk),i=1,4)
+ do i=1,4
+ bvar(i,l,j,kk)=deg2rad*bvar(i,l,j,kk)
+ enddo ! l
+ enddo ! l
+ enddo ! j
+! endif
+ enddo ! k
+
+ if (kk.lt.nbank) nbank=kk
+!d write (*,*) "read_bank ******************* kk",kk,
+!d & "nbank",nbank
+!d do kk=1,nbank
+!d print *,"kk=",kk
+!d do j=1,numch
+!d do l=2,nres-1
+!d write (*,850) (bvar(i,l,j,kk),i=1,4)
+!d enddo
+!d enddo
+!d enddo
+
+! do k=1,nbank
+! read (33,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
+! do j=1,numch
+! do l=2,nres-1
+! read (33,850) (bvar(i,l,j,k),i=1,4)
+! do i=1,4
+! bvar(i,l,j,k)=deg2rad*bvar(i,l,j,k)
+! enddo
+! enddo
+! enddo
+! enddo
+ close(icsa_bank)
+
+ 850 format (10f8.3)
+ 952 format (1x,'#',i4,' total E ',f12.3,' rmsd from N ',f8.3,i5)
+ 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
+ 902 format (1x,11x,i4,12x,1pe14.5)
+ 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2,i5)
+
+ return
+ end subroutine read_bank
+!-----------------------------------------------------------------------------
+ subroutine write_bank1(jlee)
+
+ use csa_data
+ use geometry_data, only: nres,rad2deg
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+! include 'COMMON.CHAIN'
+! include 'COMMON.GEO'
+!el local variables
+ integer :: k,i,l,j,jlee
+
+#if defined(AIX) || defined(PGI)
+ open(icsa_bank1,file=csa_bank1,position="append")
+#else
+ open(icsa_bank1,file=csa_bank1,access="append")
+#endif
+ write (icsa_bank1,900) jlee,nbank,nstep,cutdif
+ do k=1,nbank
+ write (icsa_bank1,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
+ do j=1,numch
+ do l=2,nres-1
+ write (icsa_bank1,850) (rad2deg*bvar(i,l,j,k),i=1,4)
+ enddo
+ enddo
+ enddo
+ close(icsa_bank1)
+ 850 format (10f8.3)
+ 900 format (4x,"jlee =",i5,3x,"nbank =",i5,3x,"nstep =",i10,f15.5)
+ 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
+ ' %NC ',0pf5.2,i5)
+
+ return
+ end subroutine write_bank1
+!-----------------------------------------------------------------------------
+! cartprint.f
+!-----------------------------------------------------------------------------
+! subroutine cartprint
+
+! use geometry_data, only: c
+! use energy_data, only: itype
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CHAIN'
+! include 'COMMON.INTERACT'
+! include 'COMMON.NAMES'
+! include 'COMMON.IOUNITS'
+! integer :: i
+
+! write (iout,100)
+! do i=1,nres
+! write (iout,110) restyp(itype(i)),i,c(1,i),c(2,i),&
+! c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
+! enddo
+! 100 format (//' alpha-carbon coordinates ',&
+! ' centroid coordinates'/ &
+! ' ', 6X,'X',11X,'Y',11X,'Z',&
+! 10X,'X',11X,'Y',11X,'Z')
+! 110 format (a,'(',i3,')',6f12.5)
+! return
+! end subroutine cartprint
+!-----------------------------------------------------------------------------
+! dihed_cons.F
+!-----------------------------------------------------------------------------
+ subroutine secstrp2dihc
+
+ use geometry_data
+ use energy_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.GEO'
+! include 'COMMON.BOUNDS'
+! include 'COMMON.CHAIN'
+! include 'COMMON.TORCNSTR'
+! include 'COMMON.IOUNITS'
+!el character(len=1),dimension(nres) :: secstruc !(maxres)
+!el COMMON/SECONDARYS/secstruc
+ character(len=80) :: line
+ logical :: errflag
+!el external ilen
+
+!el local variables
+ integer :: i,ii,lenpre
+
+ allocate(secstruc(nres))
+
+!dr call getenv_loc('SECPREDFIL',secpred)
+ lenpre=ilen(prefix)
+ secpred=prefix(:lenpre)//'.spred'
+
+#if defined(WINIFL) || defined(WINPGI)
+ open(isecpred,file=secpred,status='old',readonly,shared)
+#elif (defined CRAY) || (defined AIX)
+ open(isecpred,file=secpred,status='old',action='read')
+#elif (defined G77)
+ open(isecpred,file=secpred,status='old')
+#else
+ open(isecpred,file=secpred,status='old',action='read')
+#endif
+! read secondary structure prediction from JPRED here!
+! read(isecpred,'(A80)',err=100,end=100) line
+! read(line,'(f10.3)',err=110) ftors
+ read(isecpred,'(f10.3)',err=110) ftors
+
+ write (iout,*) 'FTORS factor =',ftors
+! initialize secstruc to any
+ do i=1,nres
+ secstruc(i) ='-'
+ enddo
+ ndih_constr=0
+ ndih_nconstr=0
+
+ call read_secstr_pred(isecpred,iout,errflag)
+ if (errflag) then
+ write(iout,*)'There is a problem with the list of secondary-',&
+ 'structure prediction'
+ goto 100
+ endif
+! 8/13/98 Set limits to generating the dihedral angles
+ do i=1,nres
+ phibound(1,i)=-pi
+ phibound(2,i)=pi
+ enddo
+
+ ii=0
+ do i=1,nres
+ if ( secstruc(i) .eq. 'H') then
+! Helix restraints for this residue
+ ii=ii+1
+ idih_constr(ii)=i
+ phi0(ii) = 45.0D0*deg2rad
+ drange(ii)= 5.0D0*deg2rad
+ phibound(1,i) = phi0(ii)-drange(ii)
+ phibound(2,i) = phi0(ii)+drange(ii)
+ else if (secstruc(i) .eq. 'E') then
+! strand restraints for this residue
+ ii=ii+1
+ idih_constr(ii)=i
+ phi0(ii) = 180.0D0*deg2rad
+ drange(ii)= 5.0D0*deg2rad
+ phibound(1,i) = phi0(ii)-drange(ii)
+ phibound(2,i) = phi0(ii)+drange(ii)
+ else
+! no restraints for this residue
+ ndih_nconstr=ndih_nconstr+1
+ idih_nconstr(ndih_nconstr)=i
+ endif
+ enddo
+ ndih_constr=ii
+! deallocate(secstruc)
+ return
+100 continue
+ write(iout,'(A30,A80)')'Error reading file SECPRED',secpred
+! deallocate(secstruc)
+ return
+110 continue
+ write(iout,'(A20)')'Error reading FTORS'
+! deallocate(secstruc)
+ return
+ end subroutine secstrp2dihc
+!-----------------------------------------------------------------------------
+ subroutine read_secstr_pred(jin,jout,errors)
+
+! implicit real*8 (a-h,o-z)
+! INCLUDE 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CHAIN'
+!el character(len=1),dimension(nres) :: secstruc !(maxres)
+!el COMMON/SECONDARYS/secstruc
+!el EXTERNAL ILEN
+ character(len=80) :: line,line1 !,ucase
+ logical :: errflag,errors,blankline
+
+!el local variables
+ integer :: jin,jout,iseq,ipos,ipos1,iend,il,&
+ length_of_chain
+ errors=.false.
+ read (jin,'(a)') line
+ write (jout,'(2a)') '> ',line(1:78)
+ line1=ucase(line)
+! Remember that we number full residues starting from 2, then, iseq=1 and iseq=nres
+! correspond to the end-groups. ADD to the secondary structure prediction "-" for the
+! end-groups in the input file "*.spred"
+
+ iseq=1
+ do while (index(line1,'$END').eq.0)
+! Override commented lines.
+ ipos=1
+ blankline=.false.
+ do while (.not.blankline)
+ line1=' '
+ call mykey(line,line1,ipos,blankline,errflag)
+ if (errflag) write (jout,'(2a)') &
+ 'Error when reading sequence in line: ',line
+ errors=errors .or. errflag
+ if (.not. blankline .and. .not. errflag) then
+ ipos1=2
+ iend=ilen(line1)
+!el if (iseq.le.maxres) then
+ if (line1(1:1).eq.'-' ) then
+ secstruc(iseq)=line1(1:1)
+ else if ( ( ucase(line1(1:1)).eq.'E' ) .or. &
+ ( ucase(line1(1:1)).eq.'H' ) ) then
+ secstruc(iseq)=ucase(line1(1:1))
+ else
+ errors=.true.
+ write (jout,1010) line1(1:1), iseq
+ goto 80
+ endif
+!el else
+!el errors=.true.
+!el write (jout,1000) iseq,maxres
+!el goto 80
+!el endif
+ do while (ipos1.le.iend)
+
+ iseq=iseq+1
+ il=1
+ ipos1=ipos1+1
+!el if (iseq.le.maxres) then
+ if (line1(ipos1-1:ipos1-1).eq.'-' ) then
+ secstruc(iseq)=line1(ipos1-1:ipos1-1)
+ else if((ucase(line1(ipos1-1:ipos1-1)).eq.'E').or. &
+ (ucase(line1(ipos1-1:ipos1-1)).eq.'H') ) then
+ secstruc(iseq)=ucase(line1(ipos1-1:ipos1-1))
+ else
+ errors=.true.
+ write (jout,1010) line1(ipos1-1:ipos1-1), iseq
+ goto 80
+ endif
+!el else
+!el errors=.true.
+!el write (jout,1000) iseq,maxres
+!el goto 80
+!el endif
+ enddo
+ iseq=iseq+1
+ endif
+ enddo
+ read (jin,'(a)') line
+ write (jout,'(2a)') '> ',line(1:78)
+ line1=ucase(line)
+ enddo
+
+!d write (jout,'(10a8)') (sequence(i),i=1,iseq-1)
+
+!d check whether the found length of the chain is correct.
+ length_of_chain=iseq-1
+ if (length_of_chain .ne. nres) then
+! errors=.true.
+ write (jout,'(a,i4,a,i4,a)') &
+ 'Error: the number of labels specified in $SEC_STRUC_PRED (' &
+ ,length_of_chain,') does not match with the number of residues (' &
+ ,nres,').'
+ endif
+ 80 continue
+
+ 1000 format('Error - the number of residues (',i4,&
+ ') has exceeded maximum (',i4,').')
+ 1010 format ('Error - unrecognized secondary structure label',a4,&
+ ' in position',i4)
+ return
+ end subroutine read_secstr_pred
+!#endif
+!-----------------------------------------------------------------------------
+! parmread.F
+!-----------------------------------------------------------------------------
+ subroutine parmread
+
+ use geometry_data
+ use energy_data
+ use control_data, only:maxtor,maxterm
+ use MD_data
+ use MPI_data
+!el use map_data
+ use control, only: getenv_loc
+!
+! Read the parameters of the probability distributions of the virtual-bond
+! valence angles and the side chains and energy parameters.
+!
+! Important! Energy-term weights ARE NOT read here; they are read from the
+! main input file instead, because NO defaults have yet been set for these
+! parameters.
+!
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+#ifdef MPI
+ include "mpif.h"
+ integer :: IERROR
+#endif
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CHAIN'
+! include 'COMMON.INTERACT'
+! include 'COMMON.GEO'
+! include 'COMMON.LOCAL'
+! include 'COMMON.TORSION'
+! include 'COMMON.SCCOR'
+! include 'COMMON.SCROT'
+! include 'COMMON.FFIELD'
+! include 'COMMON.NAMES'
+! include 'COMMON.SBRIDGE'
+! include 'COMMON.MD'
+! include 'COMMON.SETUP'
+ 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,LaTeX
+ real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob)
+ real(kind=8),dimension(13) :: bN
+ character(len=3) :: lancuch !,ucase
+!el local variables
+ integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm
+ integer :: maxinter,junk,kk,ii
+ real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
+ dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
+ sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
+ res1
+ integer :: ichir1,ichir2
+! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
+!el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
+!el allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
+!
+! For printing parameters after they are read set the following in the UNRES
+! C-shell script:
+!
+! setenv PRINT_PARM YES
+!
+! To print parameters in LaTeX format rather than as ASCII tables:
+!
+! setenv LATEX YES
+!
+ call getenv_loc("PRINT_PARM",lancuch)
+ lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
+ call getenv_loc("LATEX",lancuch)
+ LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
+!
+ dwa16=2.0d0**(1.0d0/6.0d0)
+ itypro=20
+! Assign virtual-bond length
+ vbl=3.8D0
+ vblinv=1.0D0/vbl
+ vblinv2=vblinv*vblinv
+!
+! 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)
+ allocate(msc(ntyp+1)) !(ntyp+1)
+ allocate(isc(ntyp+1)) !(ntyp+1)
+ allocate(restok(ntyp+1)) !(ntyp+1)
+ allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
+
+ dsc(:)=0.0d0
+ dsc_inv(:)=0.0d0
+
+#ifdef CRYST_BOND
+ read (ibond,*) vbldp0,akp,mp,ip,pstok
+ do i=1,ntyp
+ nbondterm(i)=1
+ read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(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,*) junk,vbldp0,akp,rjunk,mp,ip,pstok
+ do i=1,ntyp
+ read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
+ j=1,nbondterm(i)),msc(i),isc(i),restok(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/)')"Dynamic constants of the interaction sites:"
+ write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
+ 'inertia','Pstok'
+ write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
+ do i=1,ntyp
+ write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),&
+ vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(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)
+
+ a0thet(:)=0.0D0
+ athet(:,:,:,:)=0.0D0
+ bthet(:,:,:,:)=0.0D0
+ polthet(:,:)=0.0D0
+ gthet(:,:)=0.0D0
+ theta0(:)=0.0D0
+ sig0(:)=0.0D0
+ sigc0(:)=0.0D0
+
+#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,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
+ (bthet(j,i,1,1),j=1,2)
+ read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
+ read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
+ read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
+ sigc0(i)=sigc0(i)**2
+ enddo
+ 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
+ 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
+
+ close (ithep)
+ if (lprint) then
+ if (.not.LaTeX) 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,1,1),j=1,2),(bthet(j,i,1,1),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
+ else
+ 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
+ endif
+#else
+!
+! Read the parameters of Utheta determined from ab initio surfaces
+! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
+!
+ read (ithep,*,err=111,end=111) 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,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
+ do i=-ntyp1,-1
+ ithetyp(i)=-ithetyp(-i)
+ enddo
+
+ aa0thet(:,:,:,:)=0.0d0
+ aathet(:,:,:,:,:)=0.0d0
+ bbthet(:,:,:,:,:,:)=0.0d0
+ ccthet(:,:,:,:,:,:)=0.0d0
+ ddthet(:,:,:,:,:,:)=0.0d0
+ eethet(:,:,:,:,:,:)=0.0d0
+ ffthet(:,:,:,:,:,:,:)=0.0d0
+ ggthet(:,:,:,:,:,:,:)=0.0d0
+
+! VAR:iblock means terminally blocking group 1=non-proline 2=proline
+ do iblock=1,2
+! VAR:ntethtyp is type of theta potentials type currently 0=glycine
+! VAR:1=non-glicyne non-proline 2=proline
+! VAR:negative values for D-aminoacid
+ do i=0,nthetyp
+ do j=-nthetyp,nthetyp
+ do k=-nthetyp,nthetyp
+ read (ithep,'(6a)',end=111,err=111) res1
+ read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
+! VAR: aa0thet is variable describing the average value of Foureir
+! VAR: expansion series
+! VAR: aathet is foureir expansion in theta/2 angle for full formula
+! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
+!ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
+ read (ithep,*,end=111,err=111) &
+ (aathet(l,i,j,k,iblock),l=1,ntheterm)
+ read (ithep,*,end=111,err=111) &
+ ((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,*,end=111,err=111) &
+ (((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.
+! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
+! RECOMENTDED AFTER VERSION 3.3)
+! do i=1,nthetyp
+! do j=1,nthetyp
+! do l=1,ntheterm
+! aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
+! aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
+! enddo
+! aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
+! aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
+! enddo
+! do l=1,ntheterm
+! aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
+! enddo
+! aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
+! enddo
+! enddo
+! AND COMMENT THE LOOPS BELOW
+ 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 !iblock
+
+! TILL HERE
+! 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
+!
+ if (lprint) then
+ write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
+ do iblock=1,2
+ 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 !n
+ enddo !m
+ enddo !l
+ enddo !k
+ enddo !j
+ enddo !i
+ enddo
+ call flush(iout)
+ endif
+ write (2,*) "Start reading THETA_PDB",ithep_pdb
+ do i=1,ntyp
+! write (2,*) 'i=',i
+ read (ithep_pdb,*,err=111,end=111) &
+ a0thet(i),(athet(j,i,1,1),j=1,2),&
+ (bthet(j,i,1,1),j=1,2)
+ read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
+ read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
+ read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
+ sigc0(i)=sigc0(i)**2
+ enddo
+ 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
+ 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
+ write (2,*) "End reading THETA_PDB"
+ close (ithep_pdb)
+#endif
+ close(ithep)
+
+!-------------------------------------------
+ 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)
+
+ bsc(:,:)=0.0D0
+ nlob(:)=0
+ nlob(:)=0
+ dsc(:)=0.0D0
+ censc(:,:,:)=0.0D0
+ gaussc(:,:,:,:)=0.0D0
+
+#ifdef CRYST_SC
+!
+! Read the parameters of the probability distribution/energy expression
+! of the side chains.
+!
+ do i=1,ntyp
+ read (irotam,'(3x,i3,f8.3)',end=112,err=112) 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,*,end=112,err=112)(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,*,end=112,err=112) bsc(j,i)
+ read (irotam,*,end=112,err=112) (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
+ if (LaTeX) then
+ write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),&
+ ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
+ 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)
+ 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
+ else
+ 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,'(a)')
+ endif
+ endif
+ enddo
+ endif
+#else
+!
+! Read scrot parameters for potentials determined from all-atom AM1 calculations
+! added by Urszula Kozlowska 07/11/2007
+!
+!el Maximum number of SC local term fitting function coefficiants
+!el integer,parameter :: maxsccoef=65
+
+ allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
+
+ do i=1,ntyp
+ read (irotam,*,end=112,err=112)
+ if (i.eq.10) then
+ read (irotam,*,end=112,err=112)
+ else
+ do j=1,65
+ read(irotam,*,end=112,err=112) sc_parmin(j,i)
+ enddo
+ endif
+ enddo
+!
+! Read the parameters of the probability distribution/energy expression
+! of the side chains.
+!
+ write (2,*) "Start reading ROTAM_PDB"
+ do i=1,ntyp
+ read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) 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_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
+ ((blower(k,l,1),l=1,k),k=1,3)
+ do j=2,nlob(i)
+ read (irotam_pdb,*,end=112,err=112) bsc(j,i)
+ read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
+ ((blower(k,l,j),l=1,k),k=1,3)
+ 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
+ enddo
+ enddo
+ enddo
+ endif
+ enddo
+ close (irotam_pdb)
+ write (2,*) "End reading ROTAM_PDB"
+#endif
+ close(irotam)
+
+#ifdef CRYST_TOR
+!
+! Read torsional parameters in old format
+!
+ allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
+
+ read (itorp,*,end=113,err=113) ntortyp,nterm_old
+ if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
+ read (itorp,*,end=113,err=113) (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,*,end=113,err=113) 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,*,end=113,err=113) 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---------------------------
+ nterm(:,:,:)=0
+ nlor(:,:,:)=0
+!el---------------------------
+
+ read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
+ do i=-ntyp,-1
+ itortyp(i)=-itortyp(-i)
+ enddo
+! itortyp(ntyp1)=ntortyp
+! itortyp(-ntyp1)=-ntortyp
+ do iblock=1,2
+ write (iout,*) 'ntortyp',ntortyp
+ do i=0,ntortyp-1
+ do j=-ntortyp+1,ntortyp-1
+ read (itorp,*,end=113,err=113) 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,*,end=113,err=113) 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
+! write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
+! write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
+! v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
+ enddo
+ do k=1,nlor(i,j,iblock)
+ read (itorp,*,end=113,err=113) 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
+ write (iout,'(/a/)') 'Torsional constants:'
+ do iblock=1,2
+ do i=-ntortyp,ntortyp
+ do j=-ntortyp,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
+!elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
+!
+! 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)',end=114,err=114) 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,*,end=114,err=114) 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,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
+ ntermd_1(i,j,k,iblock))
+ read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
+ ntermd_1(i,j,k,iblock))
+ read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
+ ntermd_1(i,j,k,iblock))
+ read (itordp,*,end=114,err=114) (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,*,end=114,err=114) ((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
+! Read of Side-chain backbone correlation parameters
+! Modified 11 May 2012 by Adasko
+!CC
+!
+ read (isccor,*,end=119,err=119) nsccortyp
+
+!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)
+!-----------------------------------
+#ifdef SCCORPDB
+!el from module energy-------------
+ allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
+
+ read (isccor,*,end=119,err=119) (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
+!el from module energy---------
+ 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 l=1,maxinter
+ do i=1,nsccortyp
+ do j=1,nsccortyp
+ read (isccor,*,end=119,err=119) &
+ 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,*,end=119,err=119) 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,*,end=119,err=119) 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)
+#else
+!el from module energy-------------
+ allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
+
+ read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
+! write (iout,*) 'ntortyp',ntortyp
+ maxinter=3
+!c maxinter is maximum interaction sites
+!el from module energy---------
+ allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
+ allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
+ allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
+ allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
+!-----------------------------------
+ do l=1,maxinter
+ do i=1,nsccortyp
+ do j=1,nsccortyp
+ read (isccor,*,end=119,err=119) &
+ nterm_sccor(i,j),nlor_sccor(i,j)
+ v0ijsccor=0.0d0
+ si=-1.0d0
+
+ do k=1,nterm_sccor(i,j)
+ read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
+ v2sccor(k,l,i,j)
+ v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
+ si=-si
+ enddo
+ do k=1,nlor_sccor(i,j)
+ read (isccor,*,end=119,err=119) 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 !el ,iblock
+ enddo
+ enddo
+ enddo
+ close (isccor)
+
+#endif
+ if (lprint) then
+ write (iout,'(/a/)') 'Torsional constants:'
+ 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)
+ 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.
+!
+
+ if (lprint) then
+ write (iout,*)
+ write (iout,*) "Coefficients of the cumulants"
+ endif
+ read (ifourier,*) nloctyp
+!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)
+! el
+ b1(1,:)=0.0d0
+ b1(2,:)=0.0d0
+!--------------------------------
+
+ do i=0,nloctyp-1
+ read (ifourier,*,end=115,err=115)
+ read (ifourier,*,end=115,err=115) (bN(ii),ii=1,13)
+ if (lprint) then
+ write (iout,*) 'Type',i
+ write (iout,'(a,i2,a,f10.5)') ('bN(',ii,')=',bN(ii),ii=1,13)
+ endif
+ B1(1,i) = bN(3)
+ B1(2,i) = bN(5)
+ B1(1,-i) = bN(3)
+ B1(2,-i) = -bN(5)
+! b1(1,i)=0.0d0
+! b1(2,i)=0.0d0
+ B1tilde(1,i) = bN(3)
+ B1tilde(2,i) =-bN(5)
+ B1tilde(1,-i) =-bN(3)
+ B1tilde(2,-i) =bN(5)
+! b1tilde(1,i)=0.0d0
+! b1tilde(2,i)=0.0d0
+ B2(1,i) = bN(2)
+ B2(2,i) = bN(4)
+ B2(1,-i) =bN(2)
+ B2(2,-i) =-bN(4)
+
+! b2(1,i)=0.0d0
+! b2(2,i)=0.0d0
+ CC(1,1,i)= bN(7)
+ CC(2,2,i)=-bN(7)
+ CC(2,1,i)= bN(9)
+ CC(1,2,i)= bN(9)
+ CC(1,1,-i)= bN(7)
+ CC(2,2,-i)=-bN(7)
+ CC(2,1,-i)=-bN(9)
+ CC(1,2,-i)=-bN(9)
+! 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)=bN(7)
+ Ctilde(1,2,i)=bN(9)
+ Ctilde(2,1,i)=-bN(9)
+ Ctilde(2,2,i)=bN(7)
+ Ctilde(1,1,-i)=bN(7)
+ Ctilde(1,2,-i)=-bN(9)
+ Ctilde(2,1,-i)=bN(9)
+ Ctilde(2,2,-i)=bN(7)
+
+! 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)= bN(6)
+ DD(2,2,i)=-bN(6)
+ DD(2,1,i)= bN(8)
+ DD(1,2,i)= bN(8)
+ DD(1,1,-i)= bN(6)
+ DD(2,2,-i)=-bN(6)
+ DD(2,1,-i)=-bN(8)
+ DD(1,2,-i)=-bN(8)
+! 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)=bN(6)
+ Dtilde(1,2,i)=bN(8)
+ Dtilde(2,1,i)=-bN(8)
+ Dtilde(2,2,i)=bN(6)
+ Dtilde(1,1,-i)=bN(6)
+ Dtilde(1,2,-i)=-bN(8)
+ Dtilde(2,1,-i)=bN(8)
+ Dtilde(2,2,-i)=bN(6)
+
+! 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)= bN(10)+bN(11)
+ EE(2,2,i)=-bN(10)+bN(11)
+ EE(2,1,i)= bN(12)-bN(13)
+ EE(1,2,i)= bN(12)+bN(13)
+ EE(1,1,-i)= bN(10)+bN(11)
+ EE(2,2,-i)=-bN(10)+bN(11)
+ EE(2,1,-i)=-bN(12)+bN(13)
+ EE(1,2,-i)=-bN(12)-bN(13)
+
+! 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,*) B1(1,i),B1(2,i)
+ write (iout,*) 'B2'
+ 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,*)
+ 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,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
+ read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
+ read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
+ read (ielep,*,end=116,err=116) ((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
+! lprint=.true.
+ if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
+ ael6(i,j),ael3(i,j)
+! lprint=.false.
+ 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)
+
+ augm(:,:)=0.0D0
+ chip(:)=0.0D0
+ alp(:)=0.0D0
+ sigma0(:)=0.0D0
+ sigii(:)=0.0D0
+ rr0(:)=0.0D0
+
+!--------------------------------
+
+ read (isidep,*,end=117,err=117) ipot,expon
+ if (ipot.lt.1 .or. ipot.gt.5) then
+ write (iout,'(2a)') 'Error while reading SC interaction',&
+ 'potential file - unknown potential type.'
+#ifdef MPI
+ call MPI_Finalize(Ierror)
+#endif
+ stop
+ endif
+ expon2=expon/2
+ if(me.eq.king) &
+ 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,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
+ read (isidep,*,end=117,err=117)((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,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
+ read (isidep,*,end=117,err=117)((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,*,end=117,err=117)(eps(i,j),j=i,ntyp)
+ enddo
+ read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
+ read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
+ read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
+ read (isidep,*,end=117,err=117)(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,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
+ read (isidep,*,end=117,err=117)((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,*)"Wrong ipot"
+! 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)
+ aa(:,:)=0.0D0
+ bb(:,:)=0.0D0
+ chi(:,:)=0.0D0
+ sigma(:,:)=0.0D0
+ r0(:,:)=0.0D0
+
+!--------------------------------
+
+ 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)
+ bad(:,:)=0.0D0
+
+#ifdef OLDSCP
+!
+! Define the SC-p interaction constants (hard-coded; old style)
+!
+ do i=1,ntyp
+! "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,*,end=118,err=118) (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
+! lprint=.true.
+ if (lprint) then
+ write (iout,*) "Parameters of SC-p interactions:"
+ do i=1,ntyp
+ 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
+! lprint=.false.
+#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
+! akcm=0.0d0
+! akth=0.0d0
+! akct=0.0d0
+! v1ss=0.0d0
+! v2ss=0.0d0
+! v3ss=0.0d0
+
+ if(me.eq.king) 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
+ 111 write (iout,*) "Error reading bending energy parameters."
+ goto 999
+ 112 write (iout,*) "Error reading rotamer energy parameters."
+ goto 999
+ 113 write (iout,*) "Error reading torsional energy parameters."
+ goto 999
+ 114 write (iout,*) "Error reading double torsional energy parameters."
+ goto 999
+ 115 write (iout,*) &
+ "Error reading cumulant (multibody energy) parameters."
+ goto 999
+ 116 write (iout,*) "Error reading electrostatic energy parameters."
+ goto 999
+ 117 write (iout,*) "Error reading side chain interaction parameters."
+ goto 999
+ 118 write (iout,*) "Error reading SCp interaction parameters."
+ goto 999
+ 119 write (iout,*) "Error reading SCCOR parameters"
+ 999 continue
+#ifdef MPI
+ call MPI_Finalize(Ierror)
+#endif
+ stop
+ return
+ end subroutine parmread
+#endif
+!-----------------------------------------------------------------------------
+! printmat.f
+!-----------------------------------------------------------------------------
+ subroutine printmat(ldim,m,n,iout,key,a)
+
+ integer :: n,ldim
+ character(len=3),dimension(n) :: key
+ real(kind=8),dimension(ldim,n) :: a
+!el local variables
+ integer :: i,j,k,m,iout,nlim
+
+ do 1 i=1,n,8
+ nlim=min0(i+7,n)
+ write (iout,1000) (key(k),k=i,nlim)
+ write (iout,1020)
+ 1000 format (/5x,8(6x,a3))
+ 1020 format (/80(1h-)/)
+ do 2 j=1,n
+ write (iout,1010) key(j),(a(j,k),k=i,nlim)
+ 2 continue
+ 1 continue
+ 1010 format (a3,2x,8(f9.4))
+ return
+ end subroutine printmat
+!-----------------------------------------------------------------------------
+! readpdb.F
+!-----------------------------------------------------------------------------
+ subroutine readpdb
+! Read the PDB file and convert the peptide geometry into virtual-chain
+! geometry.
+ use geometry_data
+ use energy_data, only: itype
+ use control_data
+ use compare_data
+ use MPI_data
+ use control, only: rescode
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.LOCAL'
+! include 'COMMON.VAR'
+! include 'COMMON.CHAIN'
+! include 'COMMON.INTERACT'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.GEO'
+! include 'COMMON.NAMES'
+! include 'COMMON.CONTROL'
+! include 'COMMON.DISTFIT'
+! include 'COMMON.SETUP'
+ integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift!,ity!,&
+! ishift_pdb
+ logical :: lprn=.true.,fail
+ real(kind=8),dimension(3) :: e1,e2,e3
+ real(kind=8) :: dcj,efree_temp
+ character(len=3) :: seq,res
+ character(len=5) :: atom
+ character(len=80) :: card
+ real(kind=8),dimension(3,20) :: sccor
+ integer :: kkk,lll,icha,kupa !rescode,
+ real(kind=8) :: cou
+!el local varables
+ integer,dimension(2,maxres/3) :: hfrag_alloc
+ integer,dimension(4,maxres/3) :: bfrag_alloc
+ real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
+
+ efree_temp=0.0d0
+ ibeg=1
+ ishift1=0
+ ishift=0
+! write (2,*) "UNRES_PDB",unres_pdb
+ ires=0
+ ires_old=0
+ nres=0
+ iii=0
+ lsecondary=.false.
+ nhfrag=0
+ nbfrag=0
+!-----------------------------
+ allocate(hfrag(2,maxres/3)) !(2,maxres/3)
+ allocate(bfrag(4,maxres/3)) !(4,maxres/3)
+
+ do i=1,100000
+ read (ipdbin,'(a80)',end=10) card
+! write (iout,'(a)') card
+ if (card(:5).eq.'HELIX') then
+ nhfrag=nhfrag+1
+ lsecondary=.true.
+ read(card(22:25),*) hfrag(1,nhfrag)
+ read(card(34:37),*) hfrag(2,nhfrag)
+ endif
+ if (card(:5).eq.'SHEET') then
+ nbfrag=nbfrag+1
+ lsecondary=.true.
+ read(card(24:26),*) bfrag(1,nbfrag)
+ read(card(35:37),*) bfrag(2,nbfrag)
+!rc----------------------------------------
+!rc to be corrected !!!
+ bfrag(3,nbfrag)=bfrag(1,nbfrag)
+ bfrag(4,nbfrag)=bfrag(2,nbfrag)
+!rc----------------------------------------
+ endif
+ if (card(:3).eq.'END') then
+ goto 10
+ else if (card(:3).eq.'TER') then
+! End current chain
+ ires_old=ires+1
+ ishift1=ishift1+1
+ itype(ires_old)=ntyp1
+ ibeg=2
+! write (iout,*) "Chain ended",ires,ishift,ires_old
+ if (unres_pdb) then
+ do j=1,3
+ dc(j,ires)=sccor(j,iii)
+ enddo
+ else
+ call sccenter(ires,iii,sccor)
+ endif
+ iii=0
+ endif
+! Read free energy
+ if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
+! Fish out the ATOM cards.
+ if (index(card(1:4),'ATOM').gt.0) then
+ read (card(12:16),*) atom
+! write (iout,*) "! ",atom," !",ires
+! if (atom.eq.'CA' .or. atom.eq.'CH3') then
+ read (card(23:26),*) ires
+ read (card(18:20),'(a3)') res
+! write (iout,*) "ires",ires,ires-ishift+ishift1,
+! & " ires_old",ires_old
+! write (iout,*) "ishift",ishift," ishift1",ishift1
+! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
+ if (ires-ishift+ishift1.ne.ires_old) then
+! Calculate the CM of the preceding residue.
+! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
+ if (ibeg.eq.0) then
+! write (iout,*) "Calculating sidechain center iii",iii
+ if (unres_pdb) then
+ do j=1,3
+ dc(j,ires+nres)=sccor(j,iii)
+ enddo
+ else
+ call sccenter(ires_old,iii,sccor)
+ endif
+ iii=0
+ endif
+! Start new residue.
+ if (res.eq.'Cl-' .or. res.eq.'Na+') then
+ ires=ires_old
+ cycle
+ else if (ibeg.eq.1) then
+ write (iout,*) "BEG ires",ires
+ ishift=ires-1
+ if (res.ne.'GLY' .and. res.ne. 'ACE') then
+ ishift=ishift-1
+ itype(1)=ntyp1
+ endif
+ ires=ires-ishift+ishift1
+ ires_old=ires
+! write (iout,*) "ishift",ishift," ires",ires,&
+! " ires_old",ires_old
+ ibeg=0
+ else if (ibeg.eq.2) then
+! Start a new chain
+ ishift=-ires_old+ires-1 !!!!!
+ ishift1=ishift1-1 !!!!!
+! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
+ ires=ires-ishift+ishift1
+ ires_old=ires
+ ibeg=0
+ else
+ ishift=ishift-(ires-ishift+ishift1-ires_old-1)
+ ires=ires-ishift+ishift1
+ ires_old=ires
+ endif
+ if (res.eq.'ACE' .or. res.eq.'NHE') then
+ itype(ires)=10
+ else
+ itype(ires)=rescode(ires,res,0)
+ endif
+ else
+ ires=ires-ishift+ishift1
+ endif
+! write (iout,*) "ires_old",ires_old," ires",ires
+ if (card(27:27).eq."A" .or. card(27:27).eq."B") then
+! ishift1=ishift1+1
+ endif
+! write (2,*) "ires",ires," res ",res!," ity"!,ity
+ if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
+ res.eq.'NHE'.and.atom(:2).eq.'HN') then
+ read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
+! write (iout,*) "backbone ",atom
+#ifdef DEBUG
+ write (iout,'(2i3,2x,a,3f8.3)') &
+ ires,itype(ires),res,(c(j,ires),j=1,3)
+#endif
+ iii=iii+1
+ do j=1,3
+ sccor(j,iii)=c(j,ires)
+ enddo
+! write (*,*) card(23:27),ires,itype(ires)
+ else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
+ atom.ne.'N' .and. atom.ne.'C' .and. &
+ atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
+ atom.ne.'OXT' .and. atom(:2).ne.'3H') then
+! write (iout,*) "sidechain ",atom
+ iii=iii+1
+ read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
+ endif
+ endif
+ enddo
+ 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
+ if (ires.eq.0) return
+! Calculate dummy residue coordinates inside the "chain" of a multichain
+! system
+ nres=ires
+ do i=2,nres-1
+! write (iout,*) i,itype(i)
+ if (itype(i).eq.ntyp1) then
+! write (iout,*) "dummy",i,itype(i)
+ do j=1,3
+ c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
+! c(j,i)=(c(j,i-1)+c(j,i+1))/2
+ dc(j,i)=c(j,i)
+ enddo
+ endif
+ enddo
+! Calculate the CM of the last side chain.
+ if (iii.gt.0) then
+ if (unres_pdb) then
+ do j=1,3
+ dc(j,ires)=sccor(j,iii)
+ enddo
+ else
+ call sccenter(ires,iii,sccor)
+ endif
+ endif
+! nres=ires
+ nsup=nres
+ nstart_sup=1
+ if (itype(nres).ne.10) then
+ nres=nres+1
+ itype(nres)=ntyp1
+ if (unres_pdb) then
+! 2/15/2013 by Adam: corrected insertion of the last dummy residue
+ call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
+ if (fail) then
+ e2(1)=0.0d0
+ e2(2)=1.0d0
+ e2(3)=0.0d0
+ endif
+ do j=1,3
+ c(j,nres)=c(j,nres-1)-3.8d0*e2(j)
+ enddo
+ else
+ do j=1,3
+ dcj=c(j,nres-2)-c(j,nres-3)
+ c(j,nres)=c(j,nres-1)+dcj
+ c(j,2*nres)=c(j,nres)
+ enddo
+ endif
+ endif
+!el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
+#ifdef WHAM_RUN
+ if (nres.ne.nres0) then
+ write (iout,*) "Error: wrong parameter value: NRES=",nres,&
+ " NRES0=",nres0
+ stop "Error nres value in WHAM input"
+ endif
+#endif
+!---------------------------------
+!el reallocate tables
+! do i=1,maxres/3
+! do j=1,2
+! hfrag_alloc(j,i)=hfrag(j,i)
+! enddo
+! do j=1,4
+! bfrag_alloc(j,i)=bfrag(j,i)
+! enddo
+! enddo
+
+! deallocate(hfrag)
+! deallocate(bfrag)
+! allocate(hfrag(2,nres/3)) !(2,maxres/3)
+!el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
+!el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
+! allocate(bfrag(4,nres/3)) !(4,maxres/3)
+
+! do i=1,nhfrag
+! do j=1,2
+! hfrag(j,i)=hfrag_alloc(j,i)
+! enddo
+! enddo
+! do i=1,nbfrag
+! do j=1,4
+! bfrag(j,i)=bfrag_alloc(j,i)
+! enddo
+! enddo
+!el end reallocate tables
+!---------------------------------
+ do i=2,nres-1
+ do j=1,3
+ c(j,i+nres)=dc(j,i)
+ enddo
+ enddo
+ do j=1,3
+ c(j,nres+1)=c(j,1)
+ c(j,2*nres)=c(j,nres)
+ enddo
+ if (itype(1).eq.ntyp1) then
+ nsup=nsup-1
+ nstart_sup=2
+ if (unres_pdb) then
+! 2/15/2013 by Adam: corrected insertion of the first dummy residue
+ call refsys(2,3,4,e1,e2,e3,fail)
+ if (fail) then
+ e2(1)=0.0d0
+ e2(2)=1.0d0
+ e2(3)=0.0d0
+ endif
+ do j=1,3
+ c(j,1)=c(j,2)-3.8d0*e2(j)
+ enddo
+ else
+ do j=1,3
+ dcj=c(j,4)-c(j,3)
+ c(j,1)=c(j,2)-dcj
+ c(j,nres+1)=c(j,1)
+ enddo
+ endif
+ endif
+! Copy the coordinates to reference coordinates
+! do i=1,2*nres
+! do j=1,3
+! cref(j,i)=c(j,i)
+! enddo
+! enddo
+! Calculate internal coordinates.
+ if (lprn) then
+ write (iout,'(/a)') &
+ "Cartesian coordinates of the reference structure"
+ write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
+ "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
+ do ires=1,nres
+ write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') &
+ restyp(itype(ires)),ires,(c(j,ires),j=1,3),&
+ (c(j,ires+nres),j=1,3)
+ enddo
+ endif
+! znamy już nres wiec mozna alokowac tablice
+! Calculate internal coordinates.
+ if(me.eq.king.or..not.out1file)then
+ write (iout,'(a)') &
+ "Backbone and SC coordinates as read from the PDB"
+ do ires=1,nres
+ write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') &
+ ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),&
+ (c(j,nres+ires),j=1,3)
+ enddo
+ endif
+
+ if(.not.allocated(vbld)) then
+ allocate(vbld(2*nres))
+ do i=1,2*nres
+ vbld(i)=0.d0
+ enddo
+ endif
+ if(.not.allocated(vbld_inv)) then
+ allocate(vbld_inv(2*nres))
+ do i=1,2*nres
+ vbld_inv(i)=0.d0
+ enddo
+ endif
+!!!el
+ if(.not.allocated(theta)) then
+ allocate(theta(nres+2))
+ theta(:)=0.0d0
+ endif
+
+ if(.not.allocated(phi)) allocate(phi(nres+2))
+ if(.not.allocated(alph)) allocate(alph(nres+2))
+ if(.not.allocated(omeg)) allocate(omeg(nres+2))
+ if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
+ if(.not.allocated(phiref)) allocate(phiref(nres+2))
+ if(.not.allocated(costtab)) allocate(costtab(nres))
+ if(.not.allocated(sinttab)) allocate(sinttab(nres))
+ if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
+ if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
+ if(.not.allocated(xxref)) allocate(xxref(nres))
+ if(.not.allocated(yyref)) allocate(yyref(nres))
+ if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
+ if(.not.allocated(dc_norm)) then
+! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
+ allocate(dc_norm(3,0:2*nres+2))
+ dc_norm(:,:)=0.d0
+ endif
+
+ call int_from_cart(.true.,.false.)
+ call sc_loc_geom(.false.)
+ do i=1,nres
+ thetaref(i)=theta(i)
+ phiref(i)=phi(i)
+ enddo
+! do i=1,2*nres
+! vbld_inv(i)=0.d0
+! vbld(i)=0.d0
+! enddo
+
+ do i=1,nres-1
+ do j=1,3
+ dc(j,i)=c(j,i+1)-c(j,i)
+ dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
+ enddo
+ enddo
+ do i=2,nres-1
+ do j=1,3
+ dc(j,i+nres)=c(j,i+nres)-c(j,i)
+ dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
+ enddo
+! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
+! vbld_inv(i+nres)
+ enddo
+! call chainbuild
+! Copy the coordinates to reference coordinates
+! Splits to single chain if occurs
+
+! do i=1,2*nres
+! do j=1,3
+! cref(j,i,cou)=c(j,i)
+! enddo
+! enddo
+!
+ if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
+ if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
+ if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
+!-----------------------------
+ kkk=1
+ lll=0
+ cou=1
+ do i=1,nres
+ lll=lll+1
+!c write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
+ if (i.gt.1) then
+ if ((itype(i-1).eq.ntyp1).and.(i.gt.2)) then
+ chain_length=lll-1
+ kkk=kkk+1
+! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
+ lll=1
+ endif
+ endif
+ do j=1,3
+ cref(j,i,cou)=c(j,i)
+ cref(j,i+nres,cou)=c(j,i+nres)
+ if (i.le.nres) then
+ chain_rep(j,lll,kkk)=c(j,i)
+ chain_rep(j,lll+nres,kkk)=c(j,i+nres)
+ endif
+ enddo
+ enddo
+ write (iout,*) chain_length
+ if (chain_length.eq.0) chain_length=nres
+ do j=1,3
+ chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
+ chain_rep(j,chain_length+nres,symetr) &
+ =chain_rep(j,chain_length+nres,1)
+ enddo
+! diagnostic
+! write (iout,*) "spraw lancuchy",chain_length,symetr
+! do i=1,4
+! do kkk=1,chain_length
+! write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
+! enddo
+! enddo
+! enddiagnostic
+! makes copy of chains
+ write (iout,*) "symetr", symetr
+
+ if (symetr.gt.1) then
+ call permut(symetr)
+ nperm=1
+ do i=1,symetr
+ nperm=nperm*i
+ enddo
+ do i=1,nperm
+ write(iout,*) (tabperm(i,kkk),kkk=1,4)
+ enddo
+ do i=1,nperm
+ cou=0
+ do kkk=1,symetr
+ icha=tabperm(i,kkk)
+! write (iout,*) i,icha
+ do lll=1,chain_length
+ cou=cou+1
+ if (cou.le.nres) then
+ do j=1,3
+ kupa=mod(lll,chain_length)
+ iprzes=(kkk-1)*chain_length+lll
+ if (kupa.eq.0) kupa=chain_length
+! write (iout,*) "kupa", kupa
+ cref(j,iprzes,i)=chain_rep(j,kupa,icha)
+ cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
+ enddo
+ endif
+ enddo
+ enddo
+ enddo
+ endif
+!-koniec robienia kopii
+! diag
+ do kkk=1,nperm
+ write (iout,*) "nowa struktura", nperm
+ do i=1,nres
+ write (iout,110) restyp(itype(i)),i,cref(1,i,kkk),&
+ cref(2,i,kkk),&
+ cref(3,i,kkk),cref(1,nres+i,kkk),&
+ cref(2,nres+i,kkk),cref(3,nres+i,kkk)
+ enddo
+ 100 format (//' alpha-carbon coordinates ',&
+ ' centroid coordinates'/ &
+ ' ', 6X,'X',11X,'Y',11X,'Z', &
+ 10X,'X',11X,'Y',11X,'Z')
+ 110 format (a,'(',i3,')',6f12.5)
+
+ enddo
+!c enddiag
+ do j=1,nbfrag
+ do i=1,4
+ bfrag(i,j)=bfrag(i,j)-ishift
+ enddo
+ enddo
+
+ do j=1,nhfrag
+ do i=1,2
+ hfrag(i,j)=hfrag(i,j)-ishift
+ enddo
+ enddo
+ ishift_pdb=ishift
+
+ return
+ end subroutine readpdb
+#if !defined(WHAM_RUN) && !defined(CLUSTER)
+!-----------------------------------------------------------------------------
+! readrtns_CSA.F
+!-----------------------------------------------------------------------------
+ subroutine read_control
+!
+! Read contorl data
+!
+! use geometry_data
+ use comm_machsw
+ use energy_data
+ use control_data
+ use compare_data
+ use MCM_data
+ use map_data
+ use csa_data
+ use MD_data
+ use MPI_data
+ use random, only: random_init
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+#ifdef MP
+ use prng, only:prng_restart
+ include 'mpif.h'
+ logical :: OKRandom!, prng_restart
+ real(kind=8) :: r1
+#endif
+! include 'COMMON.IOUNITS'
+! include 'COMMON.TIME1'
+! include 'COMMON.THREAD'
+! include 'COMMON.SBRIDGE'
+! include 'COMMON.CONTROL'
+! include 'COMMON.MCM'
+! include 'COMMON.MAP'
+! include 'COMMON.HEADER'
+! include 'COMMON.CSA'
+! include 'COMMON.CHAIN'
+! include 'COMMON.MUCA'
+! include 'COMMON.MD'
+! include 'COMMON.FFIELD'
+! include 'COMMON.INTERACT'
+! include 'COMMON.SETUP'
+!el integer :: KDIAG,ICORFL,IXDR
+!el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
+ character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
+ 'EVVRSP ','Givens ','Jacobi '/),shape(diagmeth))
+! character(len=80) :: ucase
+ character(len=640) :: controlcard
+
+ real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
+
+
+ nglob_csa=0
+ eglob_csa=1d99
+ nmin_csa=0
+ read (INP,'(a)') titel
+ call card_concat(controlcard,.true.)
+! out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
+! print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
+ call reada(controlcard,'SEED',seed,0.0D0)
+ call random_init(seed)
+! Set up the time limit (caution! The time must be input in minutes!)
+ read_cart=index(controlcard,'READ_CART').gt.0
+ call readi(controlcard,'CONSTR_DIST',constr_dist,0)
+ call readi(controlcard,'SYM',symetr,1)
+ call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
+ unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
+ call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
+ call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
+ call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
+ call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
+ call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
+ call reada(controlcard,'DRMS',drms,0.1D0)
+ if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
+ write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
+ write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
+ write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
+ write (iout,'(a,f10.1)')'DRMS = ',drms
+ write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
+ write (iout,'(a,f10.1)') 'Time limit (min):',timlim
+ endif
+ call readi(controlcard,'NZ_START',nz_start,0)
+ call readi(controlcard,'NZ_END',nz_end,0)
+ call readi(controlcard,'IZ_SC',iz_sc,0)
+ timlim=60.0D0*timlim
+ safety = 60.0d0*safety
+ timem=timlim
+ modecalc=0
+ call reada(controlcard,"T_BATH",t_bath,300.0d0)
+ minim=(index(controlcard,'MINIMIZE').gt.0)
+ dccart=(index(controlcard,'CART').gt.0)
+ overlapsc=(index(controlcard,'OVERLAP').gt.0)
+ overlapsc=.not.overlapsc
+ searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
+ searchsc=.not.searchsc
+ sideadd=(index(controlcard,'SIDEADD').gt.0)
+ energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
+ outpdb=(index(controlcard,'PDBOUT').gt.0)
+ outmol2=(index(controlcard,'MOL2OUT').gt.0)
+ pdbref=(index(controlcard,'PDBREF').gt.0)
+ refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
+ indpdb=index(controlcard,'PDBSTART')
+ extconf=(index(controlcard,'EXTCONF').gt.0)
+ call readi(controlcard,'IPRINT',iprint,0)
+ call readi(controlcard,'MAXGEN',maxgen,10000)
+ call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
+ call readi(controlcard,"KDIAG",kdiag,0)
+ call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
+ if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
+ write (iout,*) "RESCALE_MODE",rescale_mode
+ split_ene=index(controlcard,'SPLIT_ENE').gt.0
+ if (index(controlcard,'REGULAR').gt.0.0D0) then
+ call reada(controlcard,'WEIDIS',weidis,0.1D0)
+ modecalc=1
+ refstr=.true.
+ endif
+ if (index(controlcard,'CHECKGRAD').gt.0) then
+ modecalc=5
+ if (index(controlcard,'CART').gt.0) then
+ icheckgrad=1
+ elseif (index(controlcard,'CARINT').gt.0) then
+ icheckgrad=2
+ else
+ icheckgrad=3
+ endif
+ elseif (index(controlcard,'THREAD').gt.0) then
+ modecalc=2
+ call readi(controlcard,'THREAD',nthread,0)
+ if (nthread.gt.0) then
+ call reada(controlcard,'WEIDIS',weidis,0.1D0)
+ else
+ if (fg_rank.eq.0) &
+ write (iout,'(a)')'A number has to follow the THREAD keyword.'
+ stop 'Error termination in Read_Control.'
+ endif
+ else if (index(controlcard,'MCMA').gt.0) then
+ modecalc=3
+ else if (index(controlcard,'MCEE').gt.0) then
+ modecalc=6
+ else if (index(controlcard,'MULTCONF').gt.0) then
+ modecalc=4
+ else if (index(controlcard,'MAP').gt.0) then
+ modecalc=7
+ call readi(controlcard,'MAP',nmap,0)
+ else if (index(controlcard,'CSA').gt.0) then
+ modecalc=8
+!rc else if (index(controlcard,'ZSCORE').gt.0) then
+!rc
+!rc ZSCORE is rm from UNRES, modecalc=9 is available
+!rc
+!rc modecalc=9
+!fcm else if (index(controlcard,'MCMF').gt.0) then
+!fmc modecalc=10
+ else if (index(controlcard,'SOFTREG').gt.0) then
+ modecalc=11
+ else if (index(controlcard,'CHECK_BOND').gt.0) then
+ modecalc=-1
+ else if (index(controlcard,'TEST').gt.0) then
+ modecalc=-2
+ else if (index(controlcard,'MD').gt.0) then
+ modecalc=12
+ else if (index(controlcard,'RE ').gt.0) then
+ modecalc=14
+ endif
+
+ lmuca=index(controlcard,'MUCA').gt.0
+ call readi(controlcard,'MUCADYN',mucadyn,0)
+ call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
+ if (lmuca .and. (me.eq.king .or. .not.out1file )) &
+ then
+ write (iout,*) 'MUCADYN=',mucadyn
+ write (iout,*) 'MUCASMOOTH=',muca_smooth
+ endif
+
+ iscode=index(controlcard,'ONE_LETTER')
+ indphi=index(controlcard,'PHI')
+ indback=index(controlcard,'BACK')
+ iranconf=index(controlcard,'RAND_CONF')
+ i2ndstr=index(controlcard,'USE_SEC_PRED')
+ gradout=index(controlcard,'GRADOUT').gt.0
+ gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
+ call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
+ if (me.eq.king .or. .not.out1file ) &
+ write (iout,*) "DISTCHAINMAX",distchainmax
+
+ if(me.eq.king.or..not.out1file) &
+ write (iout,'(2a)') diagmeth(kdiag),&
+ ' routine used to diagonalize matrices.'
+ return
+ end subroutine read_control
+!-----------------------------------------------------------------------------
+ subroutine read_REMDpar
+!
+! Read REMD settings
+!
+! use control
+! use energy
+! use geometry
+ use REMD_data
+ use MPI_data
+ use control_data, only:out1file
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.TIME1'
+! include 'COMMON.MD'
+ use MD_data
+!el #ifndef LANG0
+!el include 'COMMON.LANGEVIN'
+!el #else
+!el include 'COMMON.LANGEVIN.lang0'
+!el #endif
+! include 'COMMON.INTERACT'
+! include 'COMMON.NAMES'
+! include 'COMMON.GEO'
+! include 'COMMON.REMD'
+! include 'COMMON.CONTROL'
+! include 'COMMON.SETUP'
+! character(len=80) :: ucase
+ character(len=320) :: controlcard
+ character(len=3200) :: controlcard1
+ integer :: iremd_m_total
+!el local variables
+ integer :: i
+! real(kind=8) :: var,ene
+
+ if(me.eq.king.or..not.out1file) &
+ write (iout,*) "REMD setup"
+
+ call card_concat(controlcard,.true.)
+ call readi(controlcard,"NREP",nrep,3)
+ call readi(controlcard,"NSTEX",nstex,1000)
+ call reada(controlcard,"RETMIN",retmin,10.0d0)
+ call reada(controlcard,"RETMAX",retmax,1000.0d0)
+ mremdsync=(index(controlcard,'SYNC').gt.0)
+ call readi(controlcard,"NSYN",i_sync_step,100)
+ restart1file=(index(controlcard,'REST1FILE').gt.0)
+ traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
+ call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
+ if(max_cache_traj_use.gt.max_cache_traj) &
+ max_cache_traj_use=max_cache_traj
+ if(me.eq.king.or..not.out1file) then
+!d if (traj1file) then
+!rc caching is in testing - NTWX is not ignored
+!d write (iout,*) "NTWX value is ignored"
+!d write (iout,*) " trajectory is stored to one file by master"
+!d write (iout,*) " before exchange at NSTEX intervals"
+!d endif
+ write (iout,*) "NREP= ",nrep
+ write (iout,*) "NSTEX= ",nstex
+ write (iout,*) "SYNC= ",mremdsync
+ write (iout,*) "NSYN= ",i_sync_step
+ write (iout,*) "TRAJCACHE= ",max_cache_traj_use
+ endif
+ remd_tlist=.false.
+ allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
+ if (index(controlcard,'TLIST').gt.0) then
+ remd_tlist=.true.
+ call card_concat(controlcard1,.true.)
+ read(controlcard1,*) (remd_t(i),i=1,nrep)
+ if(me.eq.king.or..not.out1file) &
+ write (iout,*)'tlist',(remd_t(i),i=1,nrep)
+ endif
+ remd_mlist=.false.
+ if (index(controlcard,'MLIST').gt.0) then
+ remd_mlist=.true.
+ call card_concat(controlcard1,.true.)
+ read(controlcard1,*) (remd_m(i),i=1,nrep)
+ if(me.eq.king.or..not.out1file) then
+ write (iout,*)'mlist',(remd_m(i),i=1,nrep)
+ iremd_m_total=0
+ do i=1,nrep
+ iremd_m_total=iremd_m_total+remd_m(i)
+ enddo
+ write (iout,*) 'Total number of replicas ',iremd_m_total
+ endif
+ endif
+ if(me.eq.king.or..not.out1file) &
+ write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
+ return
+ end subroutine read_REMDpar
+!-----------------------------------------------------------------------------
+ subroutine read_MDpar
+!
+! Read MD settings
+!
+ use control_data, only: r_cut,rlamb,out1file
+ use energy_data
+ use geometry_data, only: pi
+ use MPI_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.TIME1'
+! include 'COMMON.MD'
+ use MD_data
+!el #ifndef LANG0
+!el include 'COMMON.LANGEVIN'
+!el #else
+!el include 'COMMON.LANGEVIN.lang0'
+!el #endif
+! include 'COMMON.INTERACT'
+! include 'COMMON.NAMES'
+! include 'COMMON.GEO'
+! include 'COMMON.SETUP'
+! include 'COMMON.CONTROL'
+! include 'COMMON.SPLITELE'
+! character(len=80) :: ucase
+ character(len=320) :: controlcard
+!el local variables
+ integer :: i
+ real(kind=8) :: eta
+
+ call card_concat(controlcard,.true.)
+ call readi(controlcard,"NSTEP",n_timestep,1000000)
+ call readi(controlcard,"NTWE",ntwe,100)
+ call readi(controlcard,"NTWX",ntwx,1000)
+ call reada(controlcard,"DT",d_time,1.0d-1)
+ call reada(controlcard,"DVMAX",dvmax,2.0d1)
+ call reada(controlcard,"DAMAX",damax,1.0d1)
+ call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
+ call readi(controlcard,"LANG",lang,0)
+ RESPA = index(controlcard,"RESPA") .gt. 0
+ call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
+ ntime_split0=ntime_split
+ call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
+ ntime_split0=ntime_split
+ call reada(controlcard,"R_CUT",r_cut,2.0d0)
+ call reada(controlcard,"LAMBDA",rlamb,0.3d0)
+ rest = index(controlcard,"REST").gt.0
+ tbf = index(controlcard,"TBF").gt.0
+ usampl = index(controlcard,"USAMPL").gt.0
+ mdpdb = index(controlcard,"MDPDB").gt.0
+ call reada(controlcard,"T_BATH",t_bath,300.0d0)
+ call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
+ call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
+ call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
+ if (count_reset_moment.eq.0) count_reset_moment=1000000000
+ call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
+ reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
+ reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
+ if (count_reset_vel.eq.0) count_reset_vel=1000000000
+ large = index(controlcard,"LARGE").gt.0
+ print_compon = index(controlcard,"PRINT_COMPON").gt.0
+ rattle = index(controlcard,"RATTLE").gt.0
+! if performing umbrella sampling, fragments constrained are read from the fragment file
+ nset=0
+ if(usampl) then
+ call read_fragments
+ endif
+
+ if(me.eq.king.or..not.out1file) then
+ write (iout,*)
+ write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
+ write (iout,*)
+ write (iout,'(a)') "The units are:"
+ write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
+ write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
+ " acceleration: angstrom/(48.9 fs)**2"
+ write (iout,'(a)') "energy: kcal/mol, temperature: K"
+ write (iout,*)
+ write (iout,'(a60,i10)') "Number of time steps:",n_timestep
+ write (iout,'(a60,f10.5,a)') &
+ "Initial time step of numerical integration:",d_time,&
+ " natural units"
+ write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
+ if (RESPA) then
+ write (iout,'(2a,i4,a)') &
+ "A-MTS algorithm used; initial time step for fast-varying",&
+ " short-range forces split into",ntime_split," steps."
+ write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
+ r_cut," lambda",rlamb
+ endif
+ write (iout,'(2a,f10.5)') &
+ "Maximum acceleration threshold to reduce the time step",&
+ "/increase split number:",damax
+ write (iout,'(2a,f10.5)') &
+ "Maximum predicted energy drift to reduce the timestep",&
+ "/increase split number:",edriftmax
+ write (iout,'(a60,f10.5)') &
+ "Maximum velocity threshold to reduce velocities:",dvmax
+ write (iout,'(a60,i10)') "Frequency of property output:",ntwe
+ write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
+ if (rattle) write (iout,'(a60)') &
+ "Rattle algorithm used to constrain the virtual bonds"
+ endif
+ reset_fricmat=1000
+ if (lang.gt.0) then
+ call reada(controlcard,"ETAWAT",etawat,0.8904d0)
+ call reada(controlcard,"RWAT",rwat,1.4d0)
+ call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
+ surfarea=index(controlcard,"SURFAREA").gt.0
+ call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
+ if(me.eq.king.or..not.out1file)then
+ write (iout,'(/a,$)') "Langevin dynamics calculation"
+ if (lang.eq.1) then
+ write (iout,'(a/)') &
+ " with direct integration of Langevin equations"
+ else if (lang.eq.2) then
+ write (iout,'(a/)') " with TINKER stochasic MD integrator"
+ else if (lang.eq.3) then
+ write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
+ else if (lang.eq.4) then
+ write (iout,'(a/)') " in overdamped mode"
+ else
+ write (iout,'(//a,i5)') &
+ "=========== ERROR: Unknown Langevin dynamics mode:",lang
+ stop
+ endif
+ write (iout,'(a60,f10.5)') "Temperature:",t_bath
+ write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
+ write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
+ write (iout,'(a60,f10.5)') &
+ "Scaling factor of the friction forces:",scal_fric
+ if (surfarea) write (iout,'(2a,i10,a)') &
+ "Friction coefficients will be scaled by solvent-accessible",&
+ " surface area every",reset_fricmat," steps."
+ endif
+! Calculate friction coefficients and bounds of stochastic forces
+ eta=6*pi*cPoise*etawat
+ if(me.eq.king.or..not.out1file) &
+ write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
+ eta
+ gamp=scal_fric*(pstok+rwat)*eta
+ stdfp=dsqrt(2*Rb*t_bath/d_time)
+ allocate(gamsc(ntyp1),stdfsc(ntyp1)) !(ntyp1)
+ do i=1,ntyp
+ gamsc(i)=scal_fric*(restok(i)+rwat)*eta
+ stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
+ enddo
+ if(me.eq.king.or..not.out1file)then
+ write (iout,'(/2a/)') &
+ "Radii of site types and friction coefficients and std's of",&
+ " stochastic forces of fully exposed sites"
+ write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
+ do i=1,ntyp
+ write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),&
+ gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
+ enddo
+ endif
+ else if (tbf) then
+ if(me.eq.king.or..not.out1file)then
+ write (iout,'(a)') "Berendsen bath calculation"
+ write (iout,'(a60,f10.5)') "Temperature:",t_bath
+ write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
+ if (reset_moment) &
+ write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
+ count_reset_moment," steps"
+ if (reset_vel) &
+ write (iout,'(a,i10,a)') &
+ "Velocities will be reset at random every",count_reset_vel,&
+ " steps"
+ endif
+ else
+ if(me.eq.king.or..not.out1file) &
+ write (iout,'(a31)') "Microcanonical mode calculation"
+ endif
+ if(me.eq.king.or..not.out1file)then
+ if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
+ if (usampl) then
+ write(iout,*) "MD running with constraints."
+ write(iout,*) "Equilibration time ", eq_time, " mtus."
+ write(iout,*) "Constraining ", nfrag," fragments."
+ write(iout,*) "Length of each fragment, weight and q0:"
+ do iset=1,nset
+ write (iout,*) "Set of restraints #",iset
+ do i=1,nfrag
+ write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
+ ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
+ enddo
+ write(iout,*) "constraints between ", npair, "fragments."
+ write(iout,*) "constraint pairs, weights and q0:"
+ do i=1,npair
+ write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
+ ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
+ enddo
+ write(iout,*) "angle constraints within ", nfrag_back,&
+ "backbone fragments."
+ write(iout,*) "fragment, weights:"
+ do i=1,nfrag_back
+ write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
+ ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
+ wfrag_back(2,i,iset),wfrag_back(3,i,iset)
+ enddo
+ enddo
+ iset=mod(kolor,nset)+1
+ endif
+ endif
+ if(me.eq.king.or..not.out1file) &
+ write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
+ return
+ end subroutine read_MDpar
+!-----------------------------------------------------------------------------
+ subroutine map_read
+
+ use map_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.MAP'
+! include 'COMMON.IOUNITS'
+ character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
+ character(len=80) :: mapcard !,ucase
+!el local variables
+ integer :: imap
+! real(kind=8) :: var,ene
+
+ do imap=1,nmap
+ read (inp,'(a)') mapcard
+ mapcard=ucase(mapcard)
+ if (index(mapcard,'PHI').gt.0) then
+ kang(imap)=1
+ else if (index(mapcard,'THE').gt.0) then
+ kang(imap)=2
+ else if (index(mapcard,'ALP').gt.0) then
+ kang(imap)=3
+ else if (index(mapcard,'OME').gt.0) then
+ kang(imap)=4
+ else
+ write(iout,'(a)')'Error - illegal variable spec in MAP card.'
+ stop 'Error - illegal variable spec in MAP card.'
+ endif
+ call readi (mapcard,'RES1',res1(imap),0)
+ call readi (mapcard,'RES2',res2(imap),0)
+ if (res1(imap).eq.0) then
+ res1(imap)=res2(imap)
+ else if (res2(imap).eq.0) then
+ res2(imap)=res1(imap)
+ endif
+ if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
+ write (iout,'(a)') &
+ 'Error - illegal definition of variable group in MAP.'
+ stop 'Error - illegal definition of variable group in MAP.'
+ endif
+ call reada(mapcard,'FROM',ang_from(imap),0.0D0)
+ call reada(mapcard,'TO',ang_to(imap),0.0D0)
+ call readi(mapcard,'NSTEP',nstep(imap),0)
+ if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
+ write (iout,'(a)') &
+ 'Illegal boundary and/or step size specification in MAP.'
+ stop 'Illegal boundary and/or step size specification in MAP.'
+ endif
+ enddo ! imap
+ return
+ end subroutine map_read
+!-----------------------------------------------------------------------------
+ subroutine csaread
+
+ use control_data, only: vdisulf
+ use csa_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.GEO'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+! include 'COMMON.CONTROL'
+! character(len=80) :: ucase
+ character(len=620) :: mcmcard
+!el local variables
+! integer :: ntf,ik,iw_pdb
+! real(kind=8) :: var,ene
+
+ call card_concat(mcmcard,.true.)
+
+ call readi(mcmcard,'NCONF',nconf,50)
+ call readi(mcmcard,'NADD',nadd,0)
+ call readi(mcmcard,'JSTART',jstart,1)
+ call readi(mcmcard,'JEND',jend,1)
+ call readi(mcmcard,'NSTMAX',nstmax,500000)
+ call readi(mcmcard,'N0',n0,1)
+ call readi(mcmcard,'N1',n1,6)
+ call readi(mcmcard,'N2',n2,4)
+ call readi(mcmcard,'N3',n3,0)
+ call readi(mcmcard,'N4',n4,0)
+ call readi(mcmcard,'N5',n5,0)
+ call readi(mcmcard,'N6',n6,10)
+ call readi(mcmcard,'N7',n7,0)
+ call readi(mcmcard,'N8',n8,0)
+ call readi(mcmcard,'N9',n9,0)
+ call readi(mcmcard,'N14',n14,0)
+ call readi(mcmcard,'N15',n15,0)
+ call readi(mcmcard,'N16',n16,0)
+ call readi(mcmcard,'N17',n17,0)
+ call readi(mcmcard,'N18',n18,0)
+
+ vdisulf=(index(mcmcard,'DYNSS').gt.0)
+
+ call readi(mcmcard,'NDIFF',ndiff,2)
+ call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
+ call readi(mcmcard,'IS1',is1,1)
+ call readi(mcmcard,'IS2',is2,8)
+ call readi(mcmcard,'NRAN0',nran0,4)
+ call readi(mcmcard,'NRAN1',nran1,2)
+ call readi(mcmcard,'IRR',irr,1)
+ call readi(mcmcard,'NSEED',nseed,20)
+ call readi(mcmcard,'NTOTAL',ntotal,10000)
+ call reada(mcmcard,'CUT1',cut1,2.0d0)
+ call reada(mcmcard,'CUT2',cut2,5.0d0)
+ call reada(mcmcard,'ESTOP',estop,-3000.0d0)
+ call readi(mcmcard,'ICMAX',icmax,3)
+ call readi(mcmcard,'IRESTART',irestart,0)
+!!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
+ ntbankm=0
+!!bankt
+ call reada(mcmcard,'DELE',dele,20.0d0)
+ call reada(mcmcard,'DIFCUT',difcut,720.0d0)
+ call readi(mcmcard,'IREF',iref,0)
+ call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
+ call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
+ call readi(mcmcard,'NCONF_IN',nconf_in,0)
+ call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
+ write (iout,*) "NCONF_IN",nconf_in
+ return
+ end subroutine csaread
+!-----------------------------------------------------------------------------
+ subroutine mcmread
+
+ use mcm_data
+ use control_data, only: MaxMoveType
+ use MD_data
+ use minim_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.MCM'
+! include 'COMMON.MCE'
+! include 'COMMON.IOUNITS'
+! character(len=80) :: ucase
+ character(len=320) :: mcmcard
+!el local variables
+ integer :: i
+! real(kind=8) :: var,ene
+
+ call card_concat(mcmcard,.true.)
+ call readi(mcmcard,'MAXACC',maxacc,100)
+ call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
+ call readi(mcmcard,'MAXTRIAL',maxtrial,100)
+ call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
+ call readi(mcmcard,'MAXREPM',maxrepm,200)
+ call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
+ call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
+ call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
+ call reada(mcmcard,'E_UP',e_up,5.0D0)
+ call reada(mcmcard,'DELTE',delte,0.1D0)
+ call readi(mcmcard,'NSWEEP',nsweep,5)
+ call readi(mcmcard,'NSTEPH',nsteph,0)
+ call readi(mcmcard,'NSTEPC',nstepc,0)
+ call reada(mcmcard,'TMIN',tmin,298.0D0)
+ call reada(mcmcard,'TMAX',tmax,298.0D0)
+ call readi(mcmcard,'NWINDOW',nwindow,0)
+ call readi(mcmcard,'PRINT_MC',print_mc,0)
+ print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
+ print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
+ ent_read=(index(mcmcard,'ENT_READ').gt.0)
+ call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
+ call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
+ call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
+ call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
+ call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
+ if (nwindow.gt.0) then
+ allocate(winstart(nwindow)) !!el (maxres)
+ allocate(winend(nwindow)) !!el
+ allocate(winlen(nwindow)) !!el
+ read (inp,*) (winstart(i),winend(i),i=1,nwindow)
+ do i=1,nwindow
+ winlen(i)=winend(i)-winstart(i)+1
+ enddo
+ endif
+ if (tmax.lt.tmin) tmax=tmin
+ if (tmax.eq.tmin) then
+ nstepc=0
+ nsteph=0
+ endif
+ if (nstepc.gt.0 .and. nsteph.gt.0) then
+ tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
+ tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
+ endif
+ allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
+! Probabilities of different move types
+ sumpro_type(0)=0.0D0
+ call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
+ call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
+ sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
+ call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
+ sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
+ call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
+ sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
+ do i=1,MaxMoveType
+ print *,'i',i,' sumprotype',sumpro_type(i)
+ sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
+ print *,'i',i,' sumprotype',sumpro_type(i)
+ enddo
+ return
+ end subroutine mcmread
+!-----------------------------------------------------------------------------
+ subroutine read_minim
+
+ use minim_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.MINIM'
+! include 'COMMON.IOUNITS'
+! character(len=80) :: ucase
+ character(len=320) :: minimcard
+!el local variables
+! integer :: ntf,ik,iw_pdb
+! real(kind=8) :: var,ene
+
+ call card_concat(minimcard,.true.)
+ call readi(minimcard,'MAXMIN',maxmin,2000)
+ call readi(minimcard,'MAXFUN',maxfun,5000)
+ call readi(minimcard,'MINMIN',minmin,maxmin)
+ call readi(minimcard,'MINFUN',minfun,maxmin)
+ call reada(minimcard,'TOLF',tolf,1.0D-2)
+ call reada(minimcard,'RTOLF',rtolf,1.0D-4)
+ print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
+ print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
+ print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
+ write (iout,'(/80(1h*)/20x,a/80(1h*))') &
+ 'Options in energy minimization:'
+ write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
+ 'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
+ 'MinMin:',MinMin,' MinFun:',MinFun,&
+ ' TolF:',TolF,' RTolF:',RTolF
+ return
+ end subroutine read_minim
+!-----------------------------------------------------------------------------
+ subroutine openunits
+
+ use energy_data, only: usampl
+ use csa_data
+ use MPI_data
+ use control_data, only:out1file
+ use control, only: getenv_loc
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+ character(len=16) :: form,nodename
+ integer :: nodelen,ierror,npos
+#endif
+! include 'COMMON.SETUP'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.MD'
+! include 'COMMON.CONTROL'
+ integer :: lenpre,lenpot,lentmp !,ilen
+!el external ilen
+ character(len=3) :: out1file_text !,ucase
+ character(len=3) :: ll
+!el external ucase
+!el local variables
+! integer :: ntf,ik,iw_pdb
+! real(kind=8) :: var,ene
+!
+! print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
+ call getenv_loc("PREFIX",prefix)
+ pref_orig = prefix
+ call getenv_loc("POT",pot)
+ call getenv_loc("DIRTMP",tmpdir)
+ call getenv_loc("CURDIR",curdir)
+ call getenv_loc("OUT1FILE",out1file_text)
+! print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
+ out1file_text=ucase(out1file_text)
+ if (out1file_text(1:1).eq."Y") then
+ out1file=.true.
+ else
+ out1file=fg_rank.gt.0
+ endif
+ lenpre=ilen(prefix)
+ lenpot=ilen(pot)
+ lentmp=ilen(tmpdir)
+ if (lentmp.gt.0) then
+ write (*,'(80(1h!))')
+ write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
+ write (*,'(80(1h!))')
+ write (*,*)"All output files will be on node /tmp directory."
+#ifdef MPI
+ call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
+ if (me.eq.king) then
+ write (*,*) "The master node is ",nodename
+ else if (fg_rank.eq.0) then
+ write (*,*) "I am the CG slave node ",nodename
+ else
+ write (*,*) "I am the FG slave node ",nodename
+ endif
+#endif
+ PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
+ lenpre = lentmp+lenpre+1
+ endif
+ entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
+! Get the names and open the input files
+#if defined(WINIFL) || defined(WINPGI)
+ open(1,file=pref_orig(:ilen(pref_orig))// &
+ '.inp',status='old',readonly,shared)
+ open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
+! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
+! Get parameter filenames and open the parameter files.
+ call getenv_loc('BONDPAR',bondname)
+ open (ibond,file=bondname,status='old',readonly,shared)
+ call getenv_loc('THETPAR',thetname)
+ open (ithep,file=thetname,status='old',readonly,shared)
+ call getenv_loc('ROTPAR',rotname)
+ open (irotam,file=rotname,status='old',readonly,shared)
+ call getenv_loc('TORPAR',torname)
+ open (itorp,file=torname,status='old',readonly,shared)
+ call getenv_loc('TORDPAR',tordname)
+ open (itordp,file=tordname,status='old',readonly,shared)
+ call getenv_loc('FOURIER',fouriername)
+ open (ifourier,file=fouriername,status='old',readonly,shared)
+ call getenv_loc('ELEPAR',elename)
+ open (ielep,file=elename,status='old',readonly,shared)
+ call getenv_loc('SIDEPAR',sidename)
+ open (isidep,file=sidename,status='old',readonly,shared)
+#elif (defined CRAY) || (defined AIX)
+ open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
+ action='read')
+! print *,"Processor",myrank," opened file 1"
+ open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
+! print *,"Processor",myrank," opened file 9"
+! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
+! Get parameter filenames and open the parameter files.
+ call getenv_loc('BONDPAR',bondname)
+ open (ibond,file=bondname,status='old',action='read')
+! print *,"Processor",myrank," opened file IBOND"
+ call getenv_loc('THETPAR',thetname)
+ open (ithep,file=thetname,status='old',action='read')
+! print *,"Processor",myrank," opened file ITHEP"
+ call getenv_loc('ROTPAR',rotname)
+ open (irotam,file=rotname,status='old',action='read')
+! print *,"Processor",myrank," opened file IROTAM"
+ call getenv_loc('TORPAR',torname)
+ open (itorp,file=torname,status='old',action='read')
+! print *,"Processor",myrank," opened file ITORP"
+ call getenv_loc('TORDPAR',tordname)
+ open (itordp,file=tordname,status='old',action='read')
+! print *,"Processor",myrank," opened file ITORDP"
+ call getenv_loc('SCCORPAR',sccorname)
+ open (isccor,file=sccorname,status='old',action='read')
+! print *,"Processor",myrank," opened file ISCCOR"
+ call getenv_loc('FOURIER',fouriername)
+ open (ifourier,file=fouriername,status='old',action='read')
+! print *,"Processor",myrank," opened file IFOURIER"
+ call getenv_loc('ELEPAR',elename)
+ open (ielep,file=elename,status='old',action='read')
+! print *,"Processor",myrank," opened file IELEP"
+ call getenv_loc('SIDEPAR',sidename)
+ open (isidep,file=sidename,status='old',action='read')
+! print *,"Processor",myrank," opened file ISIDEP"
+! print *,"Processor",myrank," opened parameter files"
+#elif (defined G77)
+ open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
+ open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
+! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
+! Get parameter filenames and open the parameter files.
+ call getenv_loc('BONDPAR',bondname)
+ open (ibond,file=bondname,status='old')
+ call getenv_loc('THETPAR',thetname)
+ open (ithep,file=thetname,status='old')
+ call getenv_loc('ROTPAR',rotname)
+ open (irotam,file=rotname,status='old')
+ call getenv_loc('TORPAR',torname)
+ open (itorp,file=torname,status='old')
+ call getenv_loc('TORDPAR',tordname)
+ open (itordp,file=tordname,status='old')
+ call getenv_loc('SCCORPAR',sccorname)
+ open (isccor,file=sccorname,status='old')
+ call getenv_loc('FOURIER',fouriername)
+ open (ifourier,file=fouriername,status='old')
+ call getenv_loc('ELEPAR',elename)
+ open (ielep,file=elename,status='old')
+ call getenv_loc('SIDEPAR',sidename)
+ open (isidep,file=sidename,status='old')
+#else
+ open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
+ readonly)
+ open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
+! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
+! Get parameter filenames and open the parameter files.
+ call getenv_loc('BONDPAR',bondname)
+ open (ibond,file=bondname,status='old',action='read')
+ call getenv_loc('THETPAR',thetname)
+ open (ithep,file=thetname,status='old',action='read')
+ call getenv_loc('ROTPAR',rotname)
+ open (irotam,file=rotname,status='old',action='read')
+ call getenv_loc('TORPAR',torname)
+ open (itorp,file=torname,status='old',action='read')
+ call getenv_loc('TORDPAR',tordname)
+ open (itordp,file=tordname,status='old',action='read')
+ call getenv_loc('SCCORPAR',sccorname)
+ open (isccor,file=sccorname,status='old',action='read')
+#ifndef CRYST_THETA
+ call getenv_loc('THETPARPDB',thetname_pdb)
+ print *,"thetname_pdb ",thetname_pdb
+ open (ithep_pdb,file=thetname_pdb,status='old',action='read')
+ print *,ithep_pdb," opened"
+#endif
+ call getenv_loc('FOURIER',fouriername)
+ open (ifourier,file=fouriername,status='old',readonly)
+ call getenv_loc('ELEPAR',elename)
+ open (ielep,file=elename,status='old',readonly)
+ call getenv_loc('SIDEPAR',sidename)
+ open (isidep,file=sidename,status='old',readonly)
+#ifndef CRYST_SC
+ call getenv_loc('ROTPARPDB',rotname_pdb)
+ open (irotam_pdb,file=rotname_pdb,status='old',action='read')
+#endif
+#endif
+#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 getenv_loc('SCPPAR',scpname)
+#if defined(WINIFL) || defined(WINPGI)
+ open (iscpp,file=scpname,status='old',readonly,shared)
+#elif (defined CRAY) || (defined AIX)
+ open (iscpp,file=scpname,status='old',action='read')
+#elif (defined G77)
+ open (iscpp,file=scpname,status='old')
+#else
+ open (iscpp,file=scpname,status='old',action='read')
+#endif
+#endif
+ call getenv_loc('PATTERN',patname)
+#if defined(WINIFL) || defined(WINPGI)
+ open (icbase,file=patname,status='old',readonly,shared)
+#elif (defined CRAY) || (defined AIX)
+ open (icbase,file=patname,status='old',action='read')
+#elif (defined G77)
+ open (icbase,file=patname,status='old')
+#else
+ open (icbase,file=patname,status='old',action='read')
+#endif
+#ifdef MPI
+! Open output file only for CG processes
+! print *,"Processor",myrank," fg_rank",fg_rank
+ if (fg_rank.eq.0) then
+
+ if (nodes.eq.1) then
+ npos=3
+ else
+ npos = dlog10(dfloat(nodes-1))+1
+ endif
+ if (npos.lt.3) npos=3
+ write (liczba,'(i1)') npos
+ form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
+ //')'
+ write (liczba,form) me
+ outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
+ liczba(:ilen(liczba))
+ intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
+ //'.int'
+ pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
+ //'.pdb'
+ mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
+ liczba(:ilen(liczba))//'.mol2'
+ statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
+ liczba(:ilen(liczba))//'.stat'
+ if (lentmp.gt.0) &
+ call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
+ //liczba(:ilen(liczba))//'.stat')
+ rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
+ //'.rst'
+ if(usampl) then
+ qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
+ liczba(:ilen(liczba))//'.const'
+ endif
+
+ endif
+#else
+ outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
+ intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
+ pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
+ mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
+ statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
+ if (lentmp.gt.0) &
+ call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
+ '.stat')
+ rest2name=prefix(:ilen(prefix))//'.rst'
+ if(usampl) then
+ qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
+ endif
+#endif
+#if defined(AIX) || defined(PGI)
+ if (me.eq.king .or. .not. out1file) &
+ open(iout,file=outname,status='unknown')
+#ifdef DEBUG
+ if (fg_rank.gt.0) then
+ write (liczba,'(i3.3)') myrank/nfgtasks
+ write (ll,'(bz,i3.3)') fg_rank
+ open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
+ status='unknown')
+ endif
+#endif
+ if(me.eq.king) then
+ open(igeom,file=intname,status='unknown',position='append')
+ open(ipdb,file=pdbname,status='unknown')
+ open(imol2,file=mol2name,status='unknown')
+ open(istat,file=statname,status='unknown',position='append')
+ else
+!1out open(iout,file=outname,status='unknown')
+ endif
+#else
+ if (me.eq.king .or. .not.out1file) &
+ open(iout,file=outname,status='unknown')
+#ifdef DEBUG
+ if (fg_rank.gt.0) then
+ write (liczba,'(i3.3)') myrank/nfgtasks
+ write (ll,'(bz,i3.3)') fg_rank
+ open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
+ status='unknown')
+ endif
+#endif
+ if(me.eq.king) then
+ open(igeom,file=intname,status='unknown',access='append')
+ open(ipdb,file=pdbname,status='unknown')
+ open(imol2,file=mol2name,status='unknown')
+ open(istat,file=statname,status='unknown',access='append')
+ else
+!1out open(iout,file=outname,status='unknown')
+ endif
+#endif
+ csa_rbank=prefix(:lenpre)//'.CSA.rbank'
+ csa_seed=prefix(:lenpre)//'.CSA.seed'
+ csa_history=prefix(:lenpre)//'.CSA.history'
+ csa_bank=prefix(:lenpre)//'.CSA.bank'
+ csa_bank1=prefix(:lenpre)//'.CSA.bank1'
+ csa_alpha=prefix(:lenpre)//'.CSA.alpha'
+ csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
+!!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
+ csa_int=prefix(:lenpre)//'.int'
+ csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
+ csa_native_int=prefix(:lenpre)//'.CSA.native.int'
+ csa_in=prefix(:lenpre)//'.CSA.in'
+! print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
+! Write file names
+ if (me.eq.king)then
+ write (iout,'(80(1h-))')
+ write (iout,'(30x,a)') "FILE ASSIGNMENT"
+ write (iout,'(80(1h-))')
+ write (iout,*) "Input file : ",&
+ pref_orig(:ilen(pref_orig))//'.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,*) "SCCOR 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))
+!el----
+#ifndef CRYST_THETA
+ write (iout,*) "Thetpdb parameter file : ",&
+ thetname_pdb(:ilen(thetname_pdb))
+#endif
+!el
+ write (iout,*) "Threading database : ",&
+ patname(:ilen(patname))
+ if (lentmp.ne.0) &
+ write (iout,*)" DIRTMP : ",&
+ tmpdir(:lentmp)
+ write (iout,'(80(1h-))')
+ endif
+ return
+ end subroutine openunits
+!-----------------------------------------------------------------------------
+ subroutine readrst
+
+ use geometry_data, only: nres,dc
+ use energy_data, only: usampl,iset
+ use MD_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CHAIN'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.MD'
+!el local variables
+ integer ::i,j
+! real(kind=8) :: var,ene
+
+ open(irest2,file=rest2name,status='unknown')
+ read(irest2,*) totT,EK,potE,totE,t_bath
+ do i=1,2*nres
+ read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
+ enddo
+ do i=1,2*nres
+ read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
+ enddo
+ if(usampl) then
+ read (irest2,*) iset
+ endif
+ close(irest2)
+ return
+ end subroutine readrst
+!-----------------------------------------------------------------------------
+ subroutine read_fragments
+
+ use energy_data
+! use geometry
+ use control_data, only:out1file
+ use MD_data
+ use MPI_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+#ifdef MPI
+ include 'mpif.h'
+#endif
+! include 'COMMON.SETUP'
+! include 'COMMON.CHAIN'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.MD'
+! include 'COMMON.CONTROL'
+!el local variables
+ integer :: i
+! real(kind=8) :: var,ene
+
+ read(inp,*) nset,nfrag,npair,nfrag_back
+
+!el from module energy
+! if(.not.allocated(mset)) allocate(mset(nset)) !(maxprocs/20)
+ if(.not.allocated(wfrag_back)) then
+ allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
+ allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
+
+ allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
+ allocate(ifrag(2,nfrag,nset)) !(2,50,maxprocs/20)
+
+ allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
+ allocate(ipair(2,npair,nset)) !(2,100,maxprocs/20)
+ endif
+
+ if(me.eq.king.or..not.out1file) &
+ write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
+ " nfrag_back",nfrag_back
+ do iset=1,nset
+ read(inp,*) mset(iset)
+ do i=1,nfrag
+ read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
+ qinfrag(i,iset)
+ if(me.eq.king.or..not.out1file) &
+ write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
+ ifrag(2,i,iset), qinfrag(i,iset)
+ enddo
+ do i=1,npair
+ read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
+ qinpair(i,iset)
+ if(me.eq.king.or..not.out1file) &
+ write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
+ ipair(2,i,iset), qinpair(i,iset)
+ enddo
+ do i=1,nfrag_back
+ read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
+ wfrag_back(3,i,iset),&
+ ifrag_back(1,i,iset),ifrag_back(2,i,iset)
+ if(me.eq.king.or..not.out1file) &
+ write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
+ wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
+ enddo
+ enddo
+ return
+ end subroutine read_fragments
+!-----------------------------------------------------------------------------
+! shift.F io_csa
+!-----------------------------------------------------------------------------
+ subroutine csa_read
+
+ use csa_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+! include 'COMMON.IOUNITS'
+!el local variables
+! integer :: ntf,ik,iw_pdb
+! real(kind=8) :: var,ene
+
+ open(icsa_in,file=csa_in,status="old",err=100)
+ read(icsa_in,*) nconf
+ read(icsa_in,*) jstart,jend
+ read(icsa_in,*) nstmax
+ read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
+ read(icsa_in,*) nran0,nran1,irr
+ read(icsa_in,*) nseed
+ read(icsa_in,*) ntotal,cut1,cut2
+ read(icsa_in,*) estop
+ read(icsa_in,*) icmax,irestart
+ read(icsa_in,*) ntbankm,dele,difcut
+ read(icsa_in,*) iref,rmscut,pnccut
+ read(icsa_in,*) ndiff
+ close(icsa_in)
+
+ return
+
+ 100 continue
+ return
+ end subroutine csa_read
+!-----------------------------------------------------------------------------
+ subroutine initial_write
+
+ use csa_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+! include 'COMMON.IOUNITS'
+!el local variables
+! integer :: ntf,ik,iw_pdb
+! real(kind=8) :: var,ene
+
+ open(icsa_seed,file=csa_seed,status="unknown")
+ write(icsa_seed,*) "seed"
+ close(31)
+#if defined(AIX) || defined(PGI)
+ open(icsa_history,file=csa_history,status="unknown",&
+ position="append")
+#else
+ open(icsa_history,file=csa_history,status="unknown",&
+ access="append")
+#endif
+ write(icsa_history,*) nconf
+ write(icsa_history,*) jstart,jend
+ write(icsa_history,*) nstmax
+ write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
+ write(icsa_history,*) nran0,nran1,irr
+ write(icsa_history,*) nseed
+ write(icsa_history,*) ntotal,cut1,cut2
+ write(icsa_history,*) estop
+ write(icsa_history,*) icmax,irestart
+ write(icsa_history,*) ntbankm,dele,difcut
+ write(icsa_history,*) iref,rmscut,pnccut
+ write(icsa_history,*) ndiff
+
+ write(icsa_history,*)
+ close(icsa_history)
+
+ open(icsa_bank1,file=csa_bank1,status="unknown")
+ write(icsa_bank1,*) 0
+ close(icsa_bank1)
+
+ return
+ end subroutine initial_write
+!-----------------------------------------------------------------------------
+ subroutine restart_write
+
+ use csa_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CSA'
+! include 'COMMON.BANK'
+!el local variables
+! integer :: ntf,ik,iw_pdb
+! real(kind=8) :: var,ene
+
+#if defined(AIX) || defined(PGI)
+ open(icsa_history,file=csa_history,position="append")
+#else
+ open(icsa_history,file=csa_history,access="append")
+#endif
+ write(icsa_history,*)
+ write(icsa_history,*) "This is restart"
+ write(icsa_history,*)
+ write(icsa_history,*) nconf
+ write(icsa_history,*) jstart,jend
+ write(icsa_history,*) nstmax
+ write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
+ write(icsa_history,*) nran0,nran1,irr
+ write(icsa_history,*) nseed
+ write(icsa_history,*) ntotal,cut1,cut2
+ write(icsa_history,*) estop
+ write(icsa_history,*) icmax,irestart
+ write(icsa_history,*) ntbankm,dele,difcut
+ write(icsa_history,*) iref,rmscut,pnccut
+ write(icsa_history,*) ndiff
+ write(icsa_history,*)
+ write(icsa_history,*) "irestart is: ", irestart
+
+ write(icsa_history,*)
+ close(icsa_history)
+
+ return
+ end subroutine restart_write
+!-----------------------------------------------------------------------------
+! test.F
+!-----------------------------------------------------------------------------
+ subroutine write_pdb(npdb,titelloc,ee)
+
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+! include 'COMMON.IOUNITS'
+ character(len=50) :: titelloc1
+ character*(*) :: titelloc
+ character(len=3) :: zahl
+ character(len=5) :: liczba5
+ real(kind=8) :: ee
+ integer :: npdb !,ilen
+!el external ilen
+!el local variables
+ integer :: lenpre
+! real(kind=8) :: var,ene
+
+ titelloc1=titelloc
+ lenpre=ilen(prefix)
+ if (npdb.lt.1000) then
+ call numstr(npdb,zahl)
+ open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
+ else
+ if (npdb.lt.10000) then
+ write(liczba5,'(i1,i4)') 0,npdb
+ else
+ write(liczba5,'(i5)') npdb
+ endif
+ open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
+ endif
+ call pdbout(ee,titelloc1,ipdb)
+ close(ipdb)
+ return
+ end subroutine write_pdb
+!-----------------------------------------------------------------------------
+! thread.F
+!-----------------------------------------------------------------------------
+ subroutine write_thread_summary
+! Thread the sequence through a database of known structures
+ use control_data, only: refstr
+! use geometry
+ use energy_data, only: n_ene_comp
+ use compare_data
+! implicit real*8 (a-h,o-z)
+! include 'DIMENSIONS'
+#ifdef MPI
+ use MPI_data !include 'COMMON.INFO'
+ include 'mpif.h'
+#endif
+! include 'COMMON.CONTROL'
+! include 'COMMON.CHAIN'
+! include 'COMMON.DBASE'
+! include 'COMMON.INTERACT'
+! include 'COMMON.VAR'
+! include 'COMMON.THREAD'
+! include 'COMMON.FFIELD'
+! include 'COMMON.SBRIDGE'
+! include 'COMMON.HEADER'
+! include 'COMMON.NAMES'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.TIME1'
+
+ integer,dimension(maxthread) :: ip
+ real(kind=8),dimension(0:n_ene) :: energia
+!el local variables
+ integer :: i,j,ii,jj,ipj,ik,kk,ist
+ real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
+
+ write (iout,'(30x,a/)') &
+ ' *********** Summary threading statistics ************'
+ write (iout,'(a)') 'Initial energies:'
+ write (iout,'(a4,2x,a12,14a14,3a8)') &
+ 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
+ 'RMSnat','NatCONT','NNCONT','RMS'
+! Energy sort patterns
+ do i=1,nthread
+ ip(i)=i
+ enddo
+ do i=1,nthread-1
+ enet=ener(n_ene-1,ip(i))
+ jj=i
+ do j=i+1,nthread
+ if (ener(n_ene-1,ip(j)).lt.enet) then
+ jj=j
+ enet=ener(n_ene-1,ip(j))
+ endif
+ enddo
+ if (jj.ne.i) then
+ ipj=ip(jj)
+ ip(jj)=ip(i)
+ ip(i)=ipj
+ endif
+ enddo
+ do ik=1,nthread
+ i=ip(ik)
+ ii=ipatt(1,i)
+ ist=nres_base(2,ii)+ipatt(2,i)
+ do kk=1,n_ene_comp
+ energia(i)=ener0(kk,i)
+ enddo
+ etot=ener0(n_ene_comp+1,i)
+ rmsnat=ener0(n_ene_comp+2,i)
+ rms=ener0(n_ene_comp+3,i)
+ frac=ener0(n_ene_comp+4,i)
+ frac_nn=ener0(n_ene_comp+5,i)
+
+ if (refstr) then
+ write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
+ i,str_nam(ii),ist+1,&
+ (energia(print_order(kk)),kk=1,nprint_ene),&
+ etot,rmsnat,frac,frac_nn,rms
+ else
+ write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
+ i,str_nam(ii),ist+1,&
+ (energia(print_order(kk)),kk=1,nprint_ene),etot
+ endif
+ enddo
+ write (iout,'(//a)') 'Final energies:'
+ write (iout,'(a4,2x,a12,17a14,3a8)') &
+ 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
+ 'RMSnat','NatCONT','NNCONT','RMS'
+ do ik=1,nthread
+ i=ip(ik)
+ ii=ipatt(1,i)
+ ist=nres_base(2,ii)+ipatt(2,i)
+ do kk=1,n_ene_comp
+ energia(kk)=ener(kk,ik)
+ enddo
+ etot=ener(n_ene_comp+1,i)
+ rmsnat=ener(n_ene_comp+2,i)
+ rms=ener(n_ene_comp+3,i)
+ frac=ener(n_ene_comp+4,i)
+ frac_nn=ener(n_ene_comp+5,i)
+ write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
+ i,str_nam(ii),ist+1,&
+ (energia(print_order(kk)),kk=1,nprint_ene),&
+ etot,rmsnat,frac,frac_nn,rms
+ enddo
+ write (iout,'(/a/)') 'IEXAM array:'
+ write (iout,'(i5)') nexcl
+ do i=1,nexcl
+ write (iout,'(2i5)') iexam(1,i),iexam(2,i)
+ enddo
+ write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
+ 'Max. time for threading step ',max_time_for_thread,&
+ 'Average time for threading step: ',ave_time_for_thread
+ return
+ end subroutine write_thread_summary
+!-----------------------------------------------------------------------------
+ subroutine write_stat_thread(ithread,ipattern,ist)
+
+ use energy_data, only: n_ene_comp
+ use compare_data
+! implicit real*8 (a-h,o-z)
+! include "DIMENSIONS"
+! include "COMMON.CONTROL"
+! include "COMMON.IOUNITS"
+! include "COMMON.THREAD"
+! include "COMMON.FFIELD"
+! include "COMMON.DBASE"
+! include "COMMON.NAMES"
+ real(kind=8),dimension(0:n_ene) :: energia
+!el local variables
+ integer :: ithread,ipattern,ist,i
+ real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
+
+#if defined(AIX) || defined(PGI)
+ open(istat,file=statname,position='append')
+#else
+ open(istat,file=statname,access='append')
+#endif
+ do i=1,n_ene_comp
+ energia(i)=ener(i,ithread)
+ enddo
+ etot=ener(n_ene_comp+1,ithread)
+ rmsnat=ener(n_ene_comp+2,ithread)
+ rms=ener(n_ene_comp+3,ithread)
+ frac=ener(n_ene_comp+4,ithread)
+ frac_nn=ener(n_ene_comp+5,ithread)
+ write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
+ ithread,str_nam(ipattern),ist+1,&
+ (energia(print_order(i)),i=1,nprint_ene),&
+ etot,rmsnat,frac,frac_nn,rms
+ close (istat)
+ return
+ end subroutine write_stat_thread
+!-----------------------------------------------------------------------------
+#endif
+!-----------------------------------------------------------------------------
+ end module io_config
projects / unres4.git / blobdiff