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