+++ /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 .not. defined(WHAM_RUN) && .not. 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) :: b
- 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) (b(ii),ii=1,13)
- if (lprint) then
- write (iout,*) 'Type',i
- write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
- endif
- B1(1,i) = b(3)
- B1(2,i) = b(5)
- B1(1,-i) = b(3)
- B1(2,-i) = -b(5)
-! b1(1,i)=0.0d0
-! b1(2,i)=0.0d0
- B1tilde(1,i) = b(3)
- B1tilde(2,i) =-b(5)
- B1tilde(1,-i) =-b(3)
- B1tilde(2,-i) =b(5)
-! b1tilde(1,i)=0.0d0
-! b1tilde(2,i)=0.0d0
- B2(1,i) = b(2)
- B2(2,i) = b(4)
- B2(1,-i) =b(2)
- B2(2,-i) =-b(4)
-
-! b2(1,i)=0.0d0
-! b2(2,i)=0.0d0
- CC(1,1,i)= b(7)
- CC(2,2,i)=-b(7)
- CC(2,1,i)= b(9)
- CC(1,2,i)= b(9)
- CC(1,1,-i)= b(7)
- CC(2,2,-i)=-b(7)
- CC(2,1,-i)=-b(9)
- CC(1,2,-i)=-b(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)=b(7)
- Ctilde(1,2,i)=b(9)
- Ctilde(2,1,i)=-b(9)
- Ctilde(2,2,i)=b(7)
- Ctilde(1,1,-i)=b(7)
- Ctilde(1,2,-i)=-b(9)
- Ctilde(2,1,-i)=b(9)
- Ctilde(2,2,-i)=b(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)= b(6)
- DD(2,2,i)=-b(6)
- DD(2,1,i)= b(8)
- DD(1,2,i)= b(8)
- DD(1,1,-i)= b(6)
- DD(2,2,-i)=-b(6)
- DD(2,1,-i)=-b(8)
- DD(1,2,-i)=-b(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)=b(6)
- Dtilde(1,2,i)=b(8)
- Dtilde(2,1,i)=-b(8)
- Dtilde(2,2,i)=b(6)
- Dtilde(1,1,-i)=b(6)
- Dtilde(1,2,-i)=-b(8)
- Dtilde(2,1,-i)=b(8)
- Dtilde(2,2,-i)=b(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)= b(10)+b(11)
- EE(2,2,i)=-b(10)+b(11)
- EE(2,1,i)= b(12)-b(13)
- EE(1,2,i)= b(12)+b(13)
- EE(1,1,-i)= b(10)+b(11)
- EE(2,2,-i)=-b(10)+b(11)
- EE(2,1,-i)=-b(12)+b(13)
- EE(1,2,-i)=-b(12)-b(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 .not. defined(WHAM_RUN) && .not. 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