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,1)),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(1) write (iout,*) 'FTORS factor =',ftors(1) ! 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 ftors(i)=ftors(1) 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:maxterm !,maxtor 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,SPLIT_FOURIERTOR real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob) real(kind=8),dimension(13) :: buse character(len=3) :: lancuch !,ucase !el local variables integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm,jj integer :: maxinter,junk,kk,ii,ncatprotparm,nkcctyp 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,epsijlip,epspeptube,epssctube,sigmapeptube, & sigmasctube integer :: ichir1,ichir2,ijunk character*3 string ! 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_nucl(ntyp_molec(2))) !(ntyp) allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp) allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp) allocate(nbondterm(ntyp)) !(ntyp) allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp) allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp) allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp) allocate(long_r_sidechain(ntyp)) allocate(short_r_sidechain(ntyp)) dsc(:)=0.0d0 dsc_inv(:)=0.0d0 #ifdef CRYST_BOND allocate(msc(ntyp+1)) !(ntyp+1) allocate(isc(ntyp+1)) !(ntyp+1) allocate(restok(ntyp+1)) !(ntyp+1) 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 mp(:)=0.0d0 ip(:)=0.0d0 msc(:,:)=0.0d0 isc(:,:)=0.0d0 allocate(msc(ntyp+1,5)) !(ntyp+1) allocate(isc(ntyp+1,5)) !(ntyp+1) allocate(restok(ntyp+1,5)) !(ntyp+1) read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp(1),ip(1),pstok(1) do i=1,ntyp_molec(1) read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),& j=1,nbondterm(i)),msc(i,1),isc(i,1),restok(i,1) 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 read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2) do i=1,ntyp_molec(2) nbondterm_nucl(i)=1 read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),isc(i,2),restok(i,2) ! dsc(i) = vbldsc0_nucl(1,i) ! if (i.eq.10) then ! dsc_inv(i)=0.0D0 ! else ! dsc_inv(i)=1.0D0/dsc(i) ! endif enddo ! read (ibond_nucl,*) junk,vbldp0_nucl,akp_nucl,rjunk,mp(2),ip(2),pstok(2) ! do i=1,ntyp_molec(2) ! read (ibond_nucl,*) nbondterm_nucl(i),(vbldsc0_nucl(j,i),& ! aksc_nucl(j,i),abond0_nucl(j,i),& ! j=1,nbondterm_nucl(i)),msc(i,2),isc(i,2),restok(i,2) ! 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 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(1),ip(1),pstok(1) do i=1,ntyp write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(i),& vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i,1),isc(i,1),restok(i,1) do j=2,nbondterm(i) write (iout,'(13x,3f10.5)') & vbldsc0(j,i),aksc(j,i),abond0(j,i) enddo enddo endif do i=1,ntyp_molec(5) read(iion,*) msc(i,5),restok(i,5) print *,msc(i,5),restok(i,5) enddo ip(5)=0.2 ! isc(5)=0.2 read (iion,*) ncatprotparm allocate(catprm(ncatprotparm,4)) do k=1,4 read (iion,*) (catprm(i,k),i=1,ncatprotparm) enddo print *, catprm ! read (iion,*) (vcatprm(k),k=1,ncatprotpram) !---------------------------------------------------- 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 allocate(liptranene(ntyp)) !C reading lipid parameters write (iout,*) "iliptranpar",iliptranpar call flush(iout) read(iliptranpar,*) pepliptran print *,pepliptran do i=1,ntyp read(iliptranpar,*) liptranene(i) print *,liptranene(i) enddo close(iliptranpar) #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,1),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,1),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,1),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,1),& 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,1),& (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,1),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 ! IF (tor_mode.eq.0) THEN 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(-nthetyp-1:nthetyp+1,& -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2)) !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2) allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,& -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2)) !(maxtheterm,-maxthetyp1:maxthetyp1,& ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2) allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,& -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2)) allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,& -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2)) allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,& -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2)) allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,& -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2)) !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,& ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2) allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,& -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2)) allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,& -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2)) !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,& ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)) read (ithep,*,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 ELSE !C here will be the apropriate recalibrating for D-aminoacid read (ithep,*,end=121,err=121) nthetyp allocate(nbend_kcc_Tb(-nthetyp:nthetyp)) allocate(v1bend_chyb(0:36,-nthetyp:nthetyp)) do i=-nthetyp+1,nthetyp-1 read (ithep,*,end=121,err=121) nbend_kcc_Tb(i) do j=0,nbend_kcc_Tb(i) read (ithep,*,end=121,err=121) ijunk,v1bend_chyb(j,i) enddo enddo if (lprint) then write (iout,'(a)') & "Parameters of the valence-only potentials" do i=-nthetyp+1,nthetyp-1 write (iout,'(2a)') "Type ",toronelet(i) do k=0,nbend_kcc_Tb(i) write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i) enddo enddo endif ENDIF ! TOR_MODE 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) !--------------- Reading theta parameters for nucleic acid------- read (ithep_nucl,*,err=111,end=111) nthetyp_nucl,ntheterm_nucl,& ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl) allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1) allocate(aa0thet_nucl(nthetyp_nucl+1,& nthetyp_nucl+1,nthetyp_nucl+1)) !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2) allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,& nthetyp_nucl+1,nthetyp_nucl+1)) !(maxtheterm,-maxthetyp1:maxthetyp1,& ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2) allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,& nthetyp_nucl+1,nthetyp_nucl+1)) allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,& nthetyp_nucl+1,nthetyp_nucl+1)) allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,& nthetyp_nucl+1,nthetyp_nucl+1)) allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,& nthetyp_nucl+1,nthetyp_nucl+1)) !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,& ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2) allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,& nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1)) allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,& nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1)) !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,& ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)) read (ithep_nucl,*,err=111,end=111) (ithetyp_nucl(i),i=1,ntyp1_molec(2)) aa0thet_nucl(:,:,:)=0.0d0 aathet_nucl(:,:,:,:)=0.0d0 bbthet_nucl(:,:,:,:,:)=0.0d0 ccthet_nucl(:,:,:,:,:)=0.0d0 ddthet_nucl(:,:,:,:,:)=0.0d0 eethet_nucl(:,:,:,:,:)=0.0d0 ffthet_nucl(:,:,:,:,:,:)=0.0d0 ggthet_nucl(:,:,:,:,:,:)=0.0d0 do i=1,nthetyp_nucl do j=1,nthetyp_nucl do k=1,nthetyp_nucl read (ithep_nucl,'(3a)',end=111,err=111) t1,t2,t3 read (ithep_nucl,*,end=111,err=111) aa0thet_nucl(i,j,k) read (ithep_nucl,*,end=111,err=111)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl) read (ithep_nucl,*,end=111,err=111) & (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), & (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), & (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), & (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl) read (ithep_nucl,*,end=111,err=111) & (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), & ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), & llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl) enddo enddo enddo !------------------------------------------- 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,1),& ' # 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 !---------reading nucleic acid parameters for rotamers------------------- allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp) do i=1,ntyp_molec(2) read (irotam_nucl,*,end=112,err=112) do j=1,9 read(irotam_nucl,*,end=112,err=112) sc_parmin_nucl(j,i) enddo enddo close(irotam_nucl) if (lprint) then write (iout,*) write (iout,*) "Base rotamer parameters" do i=1,ntyp_molec(2) write (iout,'(a)') restyp(i,2) write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9) enddo endif ! ! 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) !C !C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local !C interaction energy of the Gly, Ala, and Pro prototypes. !C read (ifourier,*) nloctyp SPLIT_FOURIERTOR = nloctyp.lt.0 nloctyp = iabs(nloctyp) !C allocate(b1(2,nres)) !(2,-maxtor:maxtor) !C allocate(b2(2,nres)) !(2,-maxtor:maxtor) !C allocate(b1tilde(2,nres)) !(2,-maxtor:maxtor) !C allocate(ctilde(2,2,nres)) !C allocate(dtilde(2,2,nres)) !(2,2,-maxtor:maxtor) !C allocate(gtb1(2,nres)) !C allocate(gtb2(2,nres)) !C allocate(cc(2,2,nres)) !C allocate(dd(2,2,nres)) !C allocate(ee(2,2,nres)) !C allocate(gtcc(2,2,nres)) !C allocate(gtdd(2,2,nres)) !C allocate(gtee(2,2,nres)) #ifdef NEWCORR allocate(itype2loc(-ntyp1:ntyp1)) allocate(iloctyp(-nloctyp:nloctyp)) allocate(bnew1(3,2,-nloctyp:nloctyp)) allocate(bnew2(3,2,-nloctyp:nloctyp)) allocate(ccnew(3,2,-nloctyp:nloctyp)) allocate(ddnew(3,2,-nloctyp:nloctyp)) allocate(e0new(3,-nloctyp:nloctyp)) allocate(eenew(2,2,2,-nloctyp:nloctyp)) allocate(bnew1tor(3,2,-nloctyp:nloctyp)) allocate(bnew2tor(3,2,-nloctyp:nloctyp)) allocate(ccnewtor(3,2,-nloctyp:nloctyp)) allocate(ddnewtor(3,2,-nloctyp:nloctyp)) allocate(e0newtor(3,-nloctyp:nloctyp)) allocate(eenewtor(2,2,2,-nloctyp:nloctyp)) read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp) read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1) itype2loc(ntyp1)=nloctyp iloctyp(nloctyp)=ntyp1 do i=1,ntyp1 itype2loc(-i)=-itype2loc(i) enddo #else allocate(iloctyp(-nloctyp:nloctyp)) allocate(itype2loc(-ntyp1:ntyp1)) iloctyp(0)=10 iloctyp(1)=9 iloctyp(2)=20 iloctyp(3)=ntyp1 #endif do i=1,nloctyp iloctyp(-i)=-iloctyp(i) enddo !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1) !c write (iout,*) "nloctyp",nloctyp, !c & " iloctyp",(iloctyp(i),i=0,nloctyp) !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1) !c write (iout,*) "nloctyp",nloctyp, !c & " iloctyp",(iloctyp(i),i=0,nloctyp) #ifdef NEWCORR do i=0,nloctyp-1 !c write (iout,*) "NEWCORR",i read (ifourier,*,end=115,err=115) do ii=1,3 do j=1,2 read (ifourier,*,end=115,err=115) bnew1(ii,j,i) enddo enddo !c write (iout,*) "NEWCORR BNEW1" !c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2) do ii=1,3 do j=1,2 read (ifourier,*,end=115,err=115) bnew2(ii,j,i) enddo enddo !c write (iout,*) "NEWCORR BNEW2" !c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2) do kk=1,3 read (ifourier,*,end=115,err=115) ccnew(kk,1,i) read (ifourier,*,end=115,err=115) ccnew(kk,2,i) enddo !c write (iout,*) "NEWCORR CCNEW" !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2) do kk=1,3 read (ifourier,*,end=115,err=115) ddnew(kk,1,i) read (ifourier,*,end=115,err=115) ddnew(kk,2,i) enddo !c write (iout,*) "NEWCORR DDNEW" !c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2) do ii=1,2 do jj=1,2 do kk=1,2 read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i) enddo enddo enddo !c write (iout,*) "NEWCORR EENEW1" !c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2) do ii=1,3 read (ifourier,*,end=115,err=115) e0new(ii,i) enddo !c write (iout,*) (e0new(ii,i),ii=1,3) enddo !c write (iout,*) "NEWCORR EENEW" do i=0,nloctyp-1 do ii=1,3 ccnew(ii,1,i)=ccnew(ii,1,i)/2 ccnew(ii,2,i)=ccnew(ii,2,i)/2 ddnew(ii,1,i)=ddnew(ii,1,i)/2 ddnew(ii,2,i)=ddnew(ii,2,i)/2 enddo enddo do i=1,nloctyp-1 do ii=1,3 bnew1(ii,1,-i)= bnew1(ii,1,i) bnew1(ii,2,-i)=-bnew1(ii,2,i) bnew2(ii,1,-i)= bnew2(ii,1,i) bnew2(ii,2,-i)=-bnew2(ii,2,i) enddo do ii=1,3 !c ccnew(ii,1,i)=ccnew(ii,1,i)/2 !c ccnew(ii,2,i)=ccnew(ii,2,i)/2 !c ddnew(ii,1,i)=ddnew(ii,1,i)/2 !c ddnew(ii,2,i)=ddnew(ii,2,i)/2 ccnew(ii,1,-i)=ccnew(ii,1,i) ccnew(ii,2,-i)=-ccnew(ii,2,i) ddnew(ii,1,-i)=ddnew(ii,1,i) ddnew(ii,2,-i)=-ddnew(ii,2,i) enddo e0new(1,-i)= e0new(1,i) e0new(2,-i)=-e0new(2,i) e0new(3,-i)=-e0new(3,i) do kk=1,2 eenew(kk,1,1,-i)= eenew(kk,1,1,i) eenew(kk,1,2,-i)=-eenew(kk,1,2,i) eenew(kk,2,1,-i)=-eenew(kk,2,1,i) eenew(kk,2,2,-i)= eenew(kk,2,2,i) enddo enddo if (lprint) then write (iout,'(a)') "Coefficients of the multibody terms" do i=-nloctyp+1,nloctyp-1 write (iout,*) "Type: ",onelet(iloctyp(i)) write (iout,*) "Coefficients of the expansion of B1" do j=1,2 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3) enddo write (iout,*) "Coefficients of the expansion of B2" do j=1,2 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3) enddo write (iout,*) "Coefficients of the expansion of C" write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3) write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3) write (iout,*) "Coefficients of the expansion of D" write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3) write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3) write (iout,*) "Coefficients of the expansion of E" write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3) do j=1,2 do k=1,2 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2) enddo enddo enddo endif IF (SPLIT_FOURIERTOR) THEN do i=0,nloctyp-1 !c write (iout,*) "NEWCORR TOR",i read (ifourier,*,end=115,err=115) do ii=1,3 do j=1,2 read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i) enddo enddo !c write (iout,*) "NEWCORR BNEW1 TOR" !c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2) do ii=1,3 do j=1,2 read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i) enddo enddo !c write (iout,*) "NEWCORR BNEW2 TOR" !c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2) do kk=1,3 read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i) read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i) enddo !c write (iout,*) "NEWCORR CCNEW TOR" !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2) do kk=1,3 read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i) read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i) enddo !c write (iout,*) "NEWCORR DDNEW TOR" !c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2) do ii=1,2 do jj=1,2 do kk=1,2 read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i) enddo enddo enddo !c write (iout,*) "NEWCORR EENEW1 TOR" !c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2) do ii=1,3 read (ifourier,*,end=115,err=115) e0newtor(ii,i) enddo !c write (iout,*) (e0newtor(ii,i),ii=1,3) enddo !c write (iout,*) "NEWCORR EENEW TOR" do i=0,nloctyp-1 do ii=1,3 ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2 ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2 ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2 ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2 enddo enddo do i=1,nloctyp-1 do ii=1,3 bnew1tor(ii,1,-i)= bnew1tor(ii,1,i) bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i) bnew2tor(ii,1,-i)= bnew2tor(ii,1,i) bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i) enddo do ii=1,3 ccnewtor(ii,1,-i)=ccnewtor(ii,1,i) ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i) ddnewtor(ii,1,-i)=ddnewtor(ii,1,i) ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i) enddo e0newtor(1,-i)= e0newtor(1,i) e0newtor(2,-i)=-e0newtor(2,i) e0newtor(3,-i)=-e0newtor(3,i) do kk=1,2 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i) eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i) eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i) eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i) enddo enddo if (lprint) then write (iout,'(a)') & "Single-body coefficients of the torsional potentials" do i=-nloctyp+1,nloctyp-1 write (iout,*) "Type: ",onelet(iloctyp(i)) write (iout,*) "Coefficients of the expansion of B1tor" do j=1,2 write (iout,'(3hB1(,i1,1h),3f10.5)') & j,(bnew1tor(k,j,i),k=1,3) enddo write (iout,*) "Coefficients of the expansion of B2tor" do j=1,2 write (iout,'(3hB2(,i1,1h),3f10.5)') & j,(bnew2tor(k,j,i),k=1,3) enddo write (iout,*) "Coefficients of the expansion of Ctor" write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3) write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3) write (iout,*) "Coefficients of the expansion of Dtor" write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3) write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3) write (iout,*) "Coefficients of the expansion of Etor" write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3) do j=1,2 do k=1,2 write (iout,'(1hE,2i1,2f10.5)') & j,k,(eenewtor(l,j,k,i),l=1,2) enddo enddo enddo endif ELSE do i=-nloctyp+1,nloctyp-1 do ii=1,3 do j=1,2 bnew1tor(ii,j,i)=bnew1(ii,j,i) enddo enddo do ii=1,3 do j=1,2 bnew2tor(ii,j,i)=bnew2(ii,j,i) enddo enddo do ii=1,3 ccnewtor(ii,1,i)=ccnew(ii,1,i) ccnewtor(ii,2,i)=ccnew(ii,2,i) ddnewtor(ii,1,i)=ddnew(ii,1,i) ddnewtor(ii,2,i)=ddnew(ii,2,i) enddo enddo ENDIF !SPLIT_FOURIER_TOR #else allocate(ccold(2,2,-nloctyp-1:nloctyp+1)) allocate(ddold(2,2,-nloctyp-1:nloctyp+1)) allocate(eeold(2,2,-nloctyp-1:nloctyp+1)) allocate(b(13,-nloctyp-1:nloctyp+1)) if (lprint) & write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)" do i=0,nloctyp-1 read (ifourier,*,end=115,err=115) read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13) if (lprint) then write (iout,*) 'Type ',onelet(iloctyp(i)) write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13) endif if (i.gt.0) then b(2,-i)= b(2,i) b(3,-i)= b(3,i) b(4,-i)=-b(4,i) b(5,-i)=-b(5,i) endif !c B1(1,i) = b(3) !c B1(2,i) = b(5) !c B1(1,-i) = b(3) !c B1(2,-i) = -b(5) !c b1(1,i)=0.0d0 !c b1(2,i)=0.0d0 !c B1tilde(1,i) = b(3) !c! B1tilde(2,i) =-b(5) !c! B1tilde(1,-i) =-b(3) !c! B1tilde(2,-i) =b(5) !c! b1tilde(1,i)=0.0d0 !c b1tilde(2,i)=0.0d0 !c B2(1,i) = b(2) !c B2(2,i) = b(4) !c B2(1,-i) =b(2) !c B2(2,-i) =-b(4) !cc B1tilde(1,i) = b(3,i) !cc B1tilde(2,i) =-b(5,i) !c B1tilde(1,-i) =-b(3,i) !c B1tilde(2,-i) =b(5,i) !cc b1tilde(1,i)=0.0d0 !cc b1tilde(2,i)=0.0d0 !cc B2(1,i) = b(2,i) !cc B2(2,i) = b(4,i) !c B2(1,-i) =b(2,i) !c B2(2,-i) =-b(4,i) !c b2(1,i)=0.0d0 !c b2(2,i)=0.0d0 CCold(1,1,i)= b(7,i) CCold(2,2,i)=-b(7,i) CCold(2,1,i)= b(9,i) CCold(1,2,i)= b(9,i) CCold(1,1,-i)= b(7,i) CCold(2,2,-i)=-b(7,i) CCold(2,1,-i)=-b(9,i) CCold(1,2,-i)=-b(9,i) !c CC(1,1,i)=0.0d0 !c CC(2,2,i)=0.0d0 !c CC(2,1,i)=0.0d0 !c CC(1,2,i)=0.0d0 !c Ctilde(1,1,i)= CCold(1,1,i) !c Ctilde(1,2,i)= CCold(1,2,i) !c Ctilde(2,1,i)=-CCold(2,1,i) !c Ctilde(2,2,i)=-CCold(2,2,i) !c CC(1,1,i)=0.0d0 !c CC(2,2,i)=0.0d0 !c CC(2,1,i)=0.0d0 !c CC(1,2,i)=0.0d0 !c Ctilde(1,1,i)= CCold(1,1,i) !c Ctilde(1,2,i)= CCold(1,2,i) !c Ctilde(2,1,i)=-CCold(2,1,i) !c Ctilde(2,2,i)=-CCold(2,2,i) !c Ctilde(1,1,i)=0.0d0 !c Ctilde(1,2,i)=0.0d0 !c Ctilde(2,1,i)=0.0d0 !c Ctilde(2,2,i)=0.0d0 DDold(1,1,i)= b(6,i) DDold(2,2,i)=-b(6,i) DDold(2,1,i)= b(8,i) DDold(1,2,i)= b(8,i) DDold(1,1,-i)= b(6,i) DDold(2,2,-i)=-b(6,i) DDold(2,1,-i)=-b(8,i) DDold(1,2,-i)=-b(8,i) !c DD(1,1,i)=0.0d0 !c DD(2,2,i)=0.0d0 !c DD(2,1,i)=0.0d0 !c DD(1,2,i)=0.0d0 !c Dtilde(1,1,i)= DD(1,1,i) !c Dtilde(1,2,i)= DD(1,2,i) !c Dtilde(2,1,i)=-DD(2,1,i) !c Dtilde(2,2,i)=-DD(2,2,i) !c Dtilde(1,1,i)=0.0d0 !c Dtilde(1,2,i)=0.0d0 !c Dtilde(2,1,i)=0.0d0 !c Dtilde(2,2,i)=0.0d0 EEold(1,1,i)= b(10,i)+b(11,i) EEold(2,2,i)=-b(10,i)+b(11,i) EEold(2,1,i)= b(12,i)-b(13,i) EEold(1,2,i)= b(12,i)+b(13,i) EEold(1,1,-i)= b(10,i)+b(11,i) EEold(2,2,-i)=-b(10,i)+b(11,i) EEold(2,1,-i)=-b(12,i)+b(13,i) EEold(1,2,-i)=-b(12,i)-b(13,i) write(iout,*) "TU DOCHODZE" print *,"JESTEM" !c ee(1,1,i)=1.0d0 !c ee(2,2,i)=1.0d0 !c ee(2,1,i)=0.0d0 !c ee(1,2,i)=0.0d0 !c ee(2,1,i)=ee(1,2,i) enddo if (lprint) then write (iout,*) write (iout,*) & "Coefficients of the cumulants (independent of valence angles)" do i=-nloctyp+1,nloctyp-1 write (iout,*) 'Type ',onelet(iloctyp(i)) write (iout,*) 'B1' write(iout,'(2f10.5)') B(3,i),B(5,i) write (iout,*) 'B2' write(iout,'(2f10.5)') B(2,i),B(4,i) write (iout,*) 'CC' do j=1,2 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i) enddo write(iout,*) 'DD' do j=1,2 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i) enddo write(iout,*) 'EE' do j=1,2 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i) enddo enddo endif #endif #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 ! IF (TOR_MODE.eq.0) THEN 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 #ifndef NEWCORR do i=1,ntyp1 itype2loc(i)=itortyp(i) enddo #endif ELSE IF (TOR_MODE.eq.1) THEN !C read valence-torsional parameters read (itorp,*,end=121,err=121) ntortyp nkcctyp=ntortyp write (iout,*) "Valence-torsional parameters read in ntortyp",& ntortyp read (itorp,*,end=121,err=121) (itortyp(i),i=1,ntyp) write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp) #ifndef NEWCORR do i=1,ntyp1 itype2loc(i)=itortyp(i) enddo #endif do i=-ntyp,-1 itortyp(i)=-itortyp(-i) enddo do i=-ntortyp+1,ntortyp-1 do j=-ntortyp+1,ntortyp-1 !C first we read the cos and sin gamma parameters read (itorp,'(13x,a)',end=121,err=121) string write (iout,*) i,j,string read (itorp,*,end=121,err=121) & nterm_kcc(j,i),nterm_kcc_Tb(j,i) !C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i) do k=1,nterm_kcc(j,i) do l=1,nterm_kcc_Tb(j,i) do ll=1,nterm_kcc_Tb(j,i) read (itorp,*,end=121,err=121) ii,jj,kk, & v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i) enddo enddo enddo enddo enddo ELSE #ifdef NEWCORR !c AL 4/8/16: Calculate coefficients from one-body parameters ntortyp=nloctyp allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1) allocate(nterm_kcc(-ntyp1:ntyp1,-ntyp1:ntyp1)) allocate(nterm_kcc_Tb(-ntyp1:ntyp1,-ntyp1:ntyp1)) allocate(v1_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1)) allocate(v2_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1)) do i=-ntyp1,ntyp1 print *,i,itortyp(i) itortyp(i)=itype2loc(i) enddo write (iout,*) & "Val-tor parameters calculated from cumulant coefficients ntortyp"& ,ntortyp do i=-ntortyp+1,ntortyp-1 do j=-ntortyp+1,ntortyp-1 nterm_kcc(j,i)=2 nterm_kcc_Tb(j,i)=3 do k=1,nterm_kcc_Tb(j,i) do l=1,nterm_kcc_Tb(j,i) v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)& +bnew1tor(k,2,i)*bnew2tor(l,2,j) v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)& +bnew1tor(k,2,i)*bnew2tor(l,1,j) enddo enddo do k=1,nterm_kcc_Tb(j,i) do l=1,nterm_kcc_Tb(j,i) #ifdef CORRCD v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j) & -ccnewtor(k,2,i)*ddnewtor(l,2,j)) v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j) & +ccnewtor(k,1,i)*ddnewtor(l,2,j)) #else v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j) & -ccnewtor(k,2,i)*ddnewtor(l,2,j)) v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j) & +ccnewtor(k,1,i)*ddnewtor(l,2,j)) #endif enddo enddo !c f(theta,gamma)=-(b21(theta)*b11(theta)+b12(theta)*b22(theta))*cos(gamma)-(b22(theta)*b11(theta)+b21(theta)*b12(theta))*sin(gamma)+(c11(theta)*d11(theta)-c12(theta)*d12(theta))*cos(2*gamma)+(c12(theta)*d11(theta)+c11(theta)*d12(theta))*sin(2*gamma) enddo enddo #else write (iout,*) "TOR_MODE>1 only with NEWCORR" stop #endif ENDIF ! TOR_MODE if (tor_mode.gt.0 .and. lprint) then !c Print valence-torsional parameters write (iout,'(a)') & "Parameters of the valence-torsional potentials" do i=-ntortyp+1,ntortyp-1 do j=-ntortyp+1,ntortyp-1 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j) write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc" do k=1,nterm_kcc(j,i) do l=1,nterm_kcc_Tb(j,i) do ll=1,nterm_kcc_Tb(j,i) write (iout,'(3i5,2f15.4)')& k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i) enddo enddo enddo enddo enddo endif #endif allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1) read (itorp_nucl,*,end=113,err=113) ntortyp_nucl ! print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2) !el from energy module--------- allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2) allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2) allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor) allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor) allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2) allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) !el--------------------------- nterm_nucl(:,:)=0 nlor_nucl(:,:)=0 !el-------------------- read (itorp_nucl,*,end=113,err=113) & (itortyp_nucl(i),i=1,ntyp_molec(2)) ! print *,itortyp_nucl(:) !c write (iout,*) 'ntortyp',ntortyp do i=1,ntortyp_nucl do j=1,ntortyp_nucl read (itorp_nucl,*,end=113,err=113) nterm_nucl(i,j),nlor_nucl(i,j) ! print *,nterm_nucl(i,j),nlor_nucl(i,j) v0ij=0.0d0 si=-1.0d0 do k=1,nterm_nucl(i,j) read (itorp_nucl,*,end=113,err=113) kk,v1_nucl(k,i,j),v2_nucl(k,i,j) v0ij=v0ij+si*v1_nucl(k,i,j) si=-si enddo do k=1,nlor_nucl(i,j) read (itorp,*,end=113,err=113) kk,vlor1_nucl(k,i,j),& vlor2_nucl(k,i,j),vlor3_nucl(k,i,j) v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2) enddo v0_nucl(i,j)=v0ij enddo enddo ! 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 l=3 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. ! ! ! 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) allocate(epslip(ntyp,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,1),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,1),sigma0(i),rr0(i),& i=1,ntyp) endif ! goto 50 !---------------------- GB or BP potential ----------------------------- case(3:4) ! 30 do i=1,ntyp ! print *,"I AM in SCELE",scelemode if (scelemode.eq.0) then 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) do i=1,ntyp read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp) enddo ! 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,1),sigma0(i),sigii(i),& chip(i),alp(i),i=1,ntyp) endif else ! print *,ntyp,"NTYP" allocate(icharge(ntyp1)) ! print *,ntyp,icharge(i) icharge(:)=0 read (isidep,*) (icharge(i),i=1,ntyp) print *,ntyp,icharge(i) ! if(.not.allocated(eps)) allocate(eps(-ntyp !c write (2,*) "icharge",(icharge(i),i=1,ntyp) allocate(alphapol(ntyp,ntyp),epshead(ntyp,ntyp),sig0head(ntyp,ntyp)) allocate(sigiso1(ntyp,ntyp),rborn(ntyp,ntyp),sigmap1(ntyp,ntyp)) allocate(sigmap2(ntyp,ntyp),sigiso2(ntyp,ntyp)) allocate(chis(ntyp,ntyp),wquad(ntyp,ntyp),chipp(ntyp,ntyp)) allocate(epsintab(ntyp,ntyp)) allocate(dtail(2,ntyp,ntyp)) allocate(alphasur(4,ntyp,ntyp),alphiso(4,ntyp,ntyp)) allocate(wqdip(2,ntyp,ntyp)) allocate(wstate(4,ntyp,ntyp)) allocate(dhead(2,2,ntyp,ntyp)) allocate(nstate(ntyp,ntyp)) allocate(debaykap(ntyp,ntyp)) if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1)) if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp) do i=1,ntyp do j=1,i ! write (*,*) "Im in ALAB", i, " ", j read(isidep,*) & eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & !6 w tej linii (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), & !6 w tej linii chis(i,j),chis(j,i), & !2 w tej linii nstate(i,j),(wstate(k,i,j),k=1,4), & !5 w tej lini - 1 integer pierwszy dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),& ! 4 w tej linii dtail(1,i,j),dtail(2,i,j), & ! 2 w tej lini epshead(i,j),sig0head(i,j), & ! 2 w tej linii rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), & ! 5 w tej linii alphapol(i,j),alphapol(j,i), & ! 2 w tej linii (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j) !8 w tej linii IF ((LaTeX).and.(i.gt.24)) then write (2,'(2a4,1h&,14(f8.2,1h&),23(f8.2,1h&))') restyp(i,1),restyp(j,1), & eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & !6 w tej linii (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), & !6 w tej linii chis(i,j),chis(j,i) !2 w tej linii endif ! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i) END DO END DO do i=1,ntyp do j=1,i IF ((LaTeX).and.(i.gt.24)) then write (2,'(2a4,1h&,14(f8.2,1h&),23(f8.2,1h&))') restyp(i,1),restyp(j,1), & dhead(1,1,i,j),dhead(2,1,i,j),& ! 2 w tej linii dtail(1,i,j),dtail(2,i,j), & ! 2 w tej lini epshead(i,j),sig0head(i,j), & ! 2 w tej linii rborn(i,j),rborn(j,i), & ! 3 w tej linii alphapol(i,j),alphapol(j,i), & ! 2 w tej linii (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j) !8 w tej linii endif END DO END DO DO i = 1, ntyp DO j = i+1, ntyp eps(i,j) = eps(j,i) sigma(i,j) = sigma(j,i) sigmap1(i,j)=sigmap1(j,i) sigmap2(i,j)=sigmap2(j,i) sigiso1(i,j)=sigiso1(j,i) sigiso2(i,j)=sigiso2(j,i) ! print *,"ATU",sigma(j,i),sigma(i,j),i,j nstate(i,j) = nstate(j,i) dtail(1,i,j) = dtail(1,j,i) dtail(2,i,j) = dtail(2,j,i) DO k = 1, 4 alphasur(k,i,j) = alphasur(k,j,i) wstate(k,i,j) = wstate(k,j,i) alphiso(k,i,j) = alphiso(k,j,i) END DO dhead(2,1,i,j) = dhead(1,1,j,i) dhead(2,2,i,j) = dhead(1,2,j,i) dhead(1,1,i,j) = dhead(2,1,j,i) dhead(1,2,i,j) = dhead(2,2,j,i) epshead(i,j) = epshead(j,i) sig0head(i,j) = sig0head(j,i) DO k = 1, 2 wqdip(k,i,j) = wqdip(k,j,i) END DO wquad(i,j) = wquad(j,i) epsintab(i,j) = epsintab(j,i) debaykap(i,j)=debaykap(j,i) ! if (epsintab(i,j).ne.1.0) print *,"WHAT?",i,j,epsintab(i,j) END DO END DO 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,1),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_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1)) if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp) allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp) if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1)) allocate(r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1) allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),& dcavtub(ntyp1)) allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),& tubetranene(ntyp1)) aa_aq(:,:)=0.0D0 bb_aq(:,:)=0.0D0 aa_lip(:,:)=0.0D0 bb_lip(:,:)=0.0D0 if (scelemode.eq.0) then chi(:,:)=0.0D0 sigma(:,:)=0.0D0 r0(:,:)=0.0D0 endif acavtub(:)=0.0d0 bcavtub(:)=0.0d0 ccavtub(:)=0.0d0 dcavtub(:)=0.0d0 sc_aa_tube_par(:)=0.0d0 sc_bb_tube_par(:)=0.0d0 !-------------------------------- do i=2,ntyp do j=1,i-1 eps(i,j)=eps(j,i) epslip(i,j)=epslip(j,i) enddo enddo if (scelemode.eq.0) then 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 endif 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) ! print *,"SIGMA ZLA?",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_aq(i,j)=epsij*rrij*rrij ! print *,"ADASKO",epsij,rrij,expon bb_aq(i,j)=-sigeps*epsij*rrij aa_aq(j,i)=aa_aq(i,j) bb_aq(j,i)=bb_aq(i,j) epsijlip=epslip(i,j) sigeps=dsign(1.0D0,epsijlip) epsijlip=dabs(epsijlip) aa_lip(i,j)=epsijlip*rrij*rrij bb_lip(i,j)=-sigeps*epsijlip*rrij aa_lip(j,i)=aa_lip(i,j) bb_lip(j,i)=bb_lip(i,j) !C write(iout,*) aa_lip if ((ipot.gt.2).and. (scelemode.eq.0)) 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,1),restyp(j,1),aa_aq(i,j),bb_aq(i,j),augm(i,j),& sigma(i,j),r0(i,j),chi(i,j),chi(j,i) endif enddo enddo allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2))) allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2)) allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2)) allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2) allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2)) allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2) allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2)) allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2) allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2) allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2) allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2) allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2) allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2) allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2) allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2)) allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2) ! augm(:,:)=0.0D0 ! chip(:)=0.0D0 ! alp(:)=0.0D0 ! sigma0(:)=0.0D0 ! sigii(:)=0.0D0 ! rr0(:)=0.0D0 read (isidep_nucl,*) ipot_nucl ! print *,"TU?!",ipot_nucl if (ipot_nucl.eq.1) then do i=1,ntyp_molec(2) do j=i,ntyp_molec(2) read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),& elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j) enddo enddo else do i=1,ntyp_molec(2) do j=i,ntyp_molec(2) read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),& chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),& elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j) enddo enddo endif ! rpp(1,1)=2**(1.0/6.0)*5.16158 do i=1,ntyp_molec(2) do j=i,ntyp_molec(2) rrij=sigma_nucl(i,j) r0_nucl(i,j)=rrij r0_nucl(j,i)=rrij rrij=rrij**expon epsij=4*eps_nucl(i,j) sigeps=dsign(1.0D0,epsij) epsij=dabs(epsij) aa_nucl(i,j)=epsij*rrij*rrij bb_nucl(i,j)=-sigeps*epsij*rrij ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij) ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- & 2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i))) enddo do j=1,i-1 aa_nucl(i,j)=aa_nucl(j,i) bb_nucl(i,j)=bb_nucl(j,i) ael3_nucl(i,j)=ael3_nucl(j,i) ael6_nucl(i,j)=ael6_nucl(j,i) ael63_nucl(i,j)=ael63_nucl(j,i) ael32_nucl(i,j)=ael32_nucl(j,i) elpp3_nucl(i,j)=elpp3_nucl(j,i) elpp6_nucl(i,j)=elpp6_nucl(j,i) elpp63_nucl(i,j)=elpp63_nucl(j,i) elpp32_nucl(i,j)=elpp32_nucl(j,i) eps_nucl(i,j)=eps_nucl(j,i) sigma_nucl(i,j)=sigma_nucl(j,i) sigmaii_nucl(i,j)=sigmaii_nucl(j,i) enddo enddo write(iout,*) "tube param" read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, & ccavtubpep,dcavtubpep,tubetranenepep sigmapeptube=sigmapeptube**6 sigeps=dsign(1.0D0,epspeptube) epspeptube=dabs(epspeptube) pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep do i=1,ntyp read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), & ccavtub(i),dcavtub(i),tubetranene(i) sigmasctube=sigmasctube**6 sigeps=dsign(1.0D0,epssctube) epssctube=dabs(epssctube) sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i) enddo !-----------------READING SC BASE POTENTIALS----------------------------- allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2))) allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2))) allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2)) allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2)) allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2))) allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2))) allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2))) allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2)) allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2))) allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2))) allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2))) allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2))) allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2))) allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2))) allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2))) allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2))) do i=1,ntyp_molec(1) do j=1,ntyp_molec(2)-1 ! without U then we will take T for U write (*,*) "Im in ", i, " ", j read(isidep_scbase,*) & eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),& chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2) write(*,*) "eps",eps_scbase(i,j) read(isidep_scbase,*) & (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), & chis_scbase(i,j,1),chis_scbase(i,j,2) read(isidep_scbase,*) & dhead_scbasei(i,j), & dhead_scbasej(i,j), & rborn_scbasei(i,j),rborn_scbasej(i,j) read(isidep_scbase,*) & (wdipdip_scbase(k,i,j),k=1,3), & (wqdip_scbase(k,i,j),k=1,2) read(isidep_scbase,*) & alphapol_scbase(i,j), & epsintab_scbase(i,j) END DO END DO allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2))) allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2))) do i=1,ntyp_molec(1) do j=1,ntyp_molec(2)-1 epsij=eps_scbase(i,j) rrij=sigma_scbase(i,j) ! r0(i,j)=rrij ! r0(j,i)=rrij rrij=rrij**expon ! epsij=eps(i,j) sigeps=dsign(1.0D0,epsij) epsij=dabs(epsij) aa_scbase(i,j)=epsij*rrij*rrij bb_scbase(i,j)=-sigeps*epsij*rrij enddo enddo !-----------------READING PEP BASE POTENTIALS------------------- allocate(eps_pepbase(ntyp_molec(2))) allocate(sigma_pepbase(ntyp_molec(2))) allocate(chi_pepbase(ntyp_molec(2),2)) allocate(chipp_pepbase(ntyp_molec(2),2)) allocate(alphasur_pepbase(4,ntyp_molec(2))) allocate(sigmap1_pepbase(ntyp_molec(2))) allocate(sigmap2_pepbase(ntyp_molec(2))) allocate(chis_pepbase(ntyp_molec(2),2)) allocate(wdipdip_pepbase(3,ntyp_molec(2))) do j=1,ntyp_molec(2)-1 ! without U then we will take T for U write (*,*) "Im in ", i, " ", j read(isidep_pepbase,*) & eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),& chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2) write(*,*) "eps",eps_pepbase(j) read(isidep_pepbase,*) & (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), & chis_pepbase(j,1),chis_pepbase(j,2) read(isidep_pepbase,*) & (wdipdip_pepbase(k,j),k=1,3) END DO allocate(aa_pepbase(ntyp_molec(2))) allocate(bb_pepbase(ntyp_molec(2))) do j=1,ntyp_molec(2)-1 epsij=eps_pepbase(j) rrij=sigma_pepbase(j) ! r0(i,j)=rrij ! r0(j,i)=rrij rrij=rrij**expon ! epsij=eps(i,j) sigeps=dsign(1.0D0,epsij) epsij=dabs(epsij) aa_pepbase(j)=epsij*rrij*rrij bb_pepbase(j)=-sigeps*epsij*rrij enddo !--------------READING SC PHOSPHATE------------------------------------- allocate(eps_scpho(ntyp_molec(1))) allocate(sigma_scpho(ntyp_molec(1))) allocate(chi_scpho(ntyp_molec(1),2)) allocate(chipp_scpho(ntyp_molec(1),2)) allocate(alphasur_scpho(4,ntyp_molec(1))) allocate(sigmap1_scpho(ntyp_molec(1))) allocate(sigmap2_scpho(ntyp_molec(1))) allocate(chis_scpho(ntyp_molec(1),2)) allocate(wqq_scpho(ntyp_molec(1))) allocate(wqdip_scpho(2,ntyp_molec(1))) allocate(alphapol_scpho(ntyp_molec(1))) allocate(epsintab_scpho(ntyp_molec(1))) allocate(dhead_scphoi(ntyp_molec(1))) allocate(rborn_scphoi(ntyp_molec(1))) allocate(rborn_scphoj(ntyp_molec(1))) allocate(alphi_scpho(ntyp_molec(1))) ! j=1 do j=1,ntyp_molec(1) ! without U then we will take T for U write (*,*) "Im in scpho ", i, " ", j read(isidep_scpho,*) & eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),& chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2) write(*,*) "eps",eps_scpho(j) read(isidep_scpho,*) & (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), & chis_scpho(j,1),chis_scpho(j,2) read(isidep_scpho,*) & (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j) read(isidep_scpho,*) & epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), & alphi_scpho(j) END DO allocate(aa_scpho(ntyp_molec(1))) allocate(bb_scpho(ntyp_molec(1))) do j=1,ntyp_molec(1) epsij=eps_scpho(j) rrij=sigma_scpho(j) ! r0(i,j)=rrij ! r0(j,i)=rrij rrij=rrij**expon ! epsij=eps(i,j) sigeps=dsign(1.0D0,epsij) epsij=dabs(epsij) aa_scpho(j)=epsij*rrij*rrij bb_scpho(j)=-sigeps*epsij*rrij enddo read(isidep_peppho,*) & eps_peppho,sigma_peppho read(isidep_peppho,*) & (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho read(isidep_peppho,*) & (wqdip_peppho(k),k=1,2) epsij=eps_peppho rrij=sigma_peppho ! r0(i,j)=rrij ! r0(j,i)=rrij rrij=rrij**expon ! epsij=eps(i,j) sigeps=dsign(1.0D0,epsij) epsij=dabs(epsij) aa_peppho=epsij*rrij*rrij bb_peppho=-sigeps*epsij*rrij 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 allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2) do i=1,ntyp_molec(2) read (iscpp_nucl,*,end=118,err=118) eps_scp_nucl(i),rscp_nucl(i) enddo do i=1,ntyp_molec(2) aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12 bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6 enddo r0pp=1.12246204830937298142*5.16158 epspp=4.95713/4 AEES=108.661 BEES=0.433246 ! ! 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 if (shield_mode.gt.0) then pi=4.0D0*datan(1.0D0) !C VSolvSphere the volume of solving sphere print *,pi,"pi" !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it !C there will be no distinction between proline peptide group and normal peptide !C group in case of shielding parameters VSolvSphere=4.0/3.0*pi*(4.50d0)**3 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3 write (iout,*) VSolvSphere,VSolvSphere_div !C long axis of side chain do i=1,ntyp long_r_sidechain(i)=vbldsc0(1,i) ! if (scelemode.eq.0) then short_r_sidechain(i)=sigma(i,i)/sqrt(2.0) if (short_r_sidechain(i).eq.0.0) short_r_sidechain(i)=0.2 ! else ! short_r_sidechain(i)=sigma(i,i) ! endif write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),& sigma0(i) enddo buff_shield=1.0d0 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" go to 999 121 write (iout,*) "Error in Czybyshev 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,istype use control_data use compare_data use MPI_data ! use control, only: rescode,sugarcode ! 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,k!,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,res2 character(len=5) :: atom character(len=80) :: card real(kind=8),dimension(3,20) :: sccor integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode, integer :: isugar,molecprev,firstion character*1 :: sugar real(kind=8) :: cou real(kind=8),dimension(3) :: ccc !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) real(kind=8),dimension(:,:), allocatable :: c_temporary integer,dimension(:,:) , allocatable :: itype_temporary integer,dimension(:),allocatable :: istype_temp efree_temp=0.0d0 ibeg=1 ishift1=0 ishift=0 molecule=1 counter=0 ! write (2,*) "UNRES_PDB",unres_pdb ires=0 ires_old=0 #ifdef WHAM_RUN do i=1,nres do j=1,5 itype(i,j)=0 enddo enddo #endif nres=0 iii=0 lsecondary=.false. nhfrag=0 nbfrag=0 do j=1,5 nres_molec(j)=0 enddo !----------------------------- allocate(hfrag(2,maxres/3)) !(2,maxres/3) allocate(bfrag(4,maxres/3)) !(4,maxres/3) if(.not. allocated(istype)) allocate(istype(maxres)) 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+2 ishift=ishift+1 ishift1=ishift1+1 itype(ires_old,molecule)=ntyp1_molec(molecule) itype(ires_old-1,molecule)=ntyp1_molec(molecule) nres_molec(molecule)=nres_molec(molecule)+2 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) ! iii=0 endif iii=0 endif ! Read free energy if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp ! Fish out the ATOM cards. ! write(iout,*) 'card',card(1:20) ! print *,"ATU ",card(1:6), CARD(3:6) ! print *,card 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+ishift1-ishift-1)=sccor(j,iii) enddo else call sccenter(ires_old,iii,sccor) endif !unres_pdb iii=0 endif !ind_pdb ! 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,1)=ntyp1 nres_molec(molecule)=nres_molec(molecule)+1 endif ! Gly 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 ! print *,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 ! Na Cl ! print *,'atom',ires,atom if (res.eq.'ACE' .or. res.eq.'NHE') then itype(ires,1)=10 else if (atom.eq.'CA '.or.atom.eq.'N ') then molecule=1 itype(ires,molecule)=rescode(ires,res,0,molecule) firstion=0 ! nres_molec(molecule)=nres_molec(molecule)+1 else molecule=2 res2=res(2:3) itype(ires,molecule)=rescode(ires,res2,0,molecule) ! nres_molec(molecule)=nres_molec(molecule)+1 read (card(19:19),'(a1)') sugar isugar=sugarcode(sugar,ires) ! if (ibeg.eq.1) then ! istype(1)=isugar ! else istype(ires)=isugar ! print *,"ires=",ires,istype(ires) ! endif endif ! atom.eq.CA endif !ACE else ires=ires-ishift+ishift1 endif !ires_old ! 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) ! print *,ires,ishift,ishift1 ! write (iout,*) "backbone ",atom #ifdef DEBUG write (iout,'(2i3,2x,a,3f8.3)') & ires,itype(ires,1),res,(c(j,ires),j=1,3) #endif iii=iii+1 nres_molec(molecule)=nres_molec(molecule)+1 do j=1,3 sccor(j,iii)=c(j,ires) enddo else if (.not.unres_pdb .and. (atom.eq."C1'" .or. & atom.eq."C2'" .or. atom.eq."C3'" & .or. atom.eq."C4'" .or. atom.eq."O4'")) then read(card(31:54),'(3f8.3)') (ccc(j),j=1,3) !c write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3) ! print *,ires,ishift,ishift1 counter=counter+1 ! iii=iii+1 ! do j=1,3 ! sccor(j,iii)=c(j,ires) ! enddo do j=1,3 c(j,ires)=c(j,ires)+ccc(j)/5.0 enddo print *,counter,molecule if (counter.eq.5) then ! iii=iii+1 nres_molec(molecule)=nres_molec(molecule)+1 firstion=0 ! do j=1,3 ! sccor(j,iii)=c(j,ires) ! enddo counter=0 endif ! print *, "ATOM",atom(1:3) else if (atom.eq."C5'") then read (card(19:19),'(a1)') sugar isugar=sugarcode(sugar,ires) if (ibeg.eq.1) then istype(1)=isugar else istype(ires)=isugar ! print *,ires,istype(ires) endif if (unres_pdb) then molecule=2 ! print *,"nres_molec(molecule)",nres_molec(molecule),ires read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3) nres_molec(molecule)=nres_molec(molecule)+1 print *,"nres_molec(molecule)",nres_molec(molecule),ires else iii=iii+1 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3) endif else if ((atom.eq."C1'").and.unres_pdb) then iii=iii+1 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3) ! write (*,*) card(23:27),ires,itype(ires,1) 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' & .and. atom.ne.'P '.and. & atom(1:1).ne.'H' .and. & atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'& ) then ! write (iout,*) "sidechain ",atom ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3) if ((molecule.ne.2).or.(atom(3:3).ne."'")) then ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3) iii=iii+1 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3) endif endif ! print *,"IONS",ions,card(1:6) else if ((ions).and.(card(1:6).eq.'HETATM')) then if (firstion.eq.0) then firstion=1 if (unres_pdb) then do j=1,3 dc(j,ires)=sccor(j,iii) enddo else call sccenter(ires,iii,sccor) endif ! unres_pdb endif !firstion read (card(12:16),*) atom ! print *,"HETATOM", atom read (card(18:20),'(a3)') res if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.& (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG') & .or.(atom(1:2).eq.'K ')) & then ires=ires+1 if (molecule.ne.5) molecprev=molecule molecule=5 nres_molec(molecule)=nres_molec(molecule)+1 print *,"HERE",nres_molec(molecule) res=res(2:3)//' ' itype(ires,molecule)=rescode(ires,res,0,molecule) read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3) endif! NA endif !atom 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 if (((ires_old.ne.ires).and.(molecule.ne.5)) & ) & nres_molec(molecule)=nres_molec(molecule)-2 print *,'I have',nres, nres_molec(:) do k=1,4 ! ions are without dummy if (nres_molec(k).eq.0) cycle do i=2,nres-1 ! write (iout,*) i,itype(i,1) ! if (itype(i,1).eq.ntyp1) then ! write (iout,*) "dummy",i,itype(i,1) ! 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 if (itype(i,k).eq.ntyp1_molec(k)) then if (itype(i+1,k).eq.ntyp1_molec(k)) then if (itype(i+2,k).eq.0) then ! print *,"masz sieczke" do j=1,5 if (itype(i+2,j).ne.0) then itype(i+1,k)=0 itype(i+1,j)=ntyp1_molec(j) nres_molec(k)=nres_molec(k)-1 nres_molec(j)=nres_molec(j)+1 go to 3331 endif enddo 3331 continue endif ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false ! if (unres_pdb) then ! 2/15/2013 by Adam: corrected insertion of the last dummy residue ! print *,i,'tu dochodze' ! call refsys(i-3,i-2,i-1,e1,e2,e3,fail) ! if (fail) then ! e2(1)=0.0d0 ! e2(2)=1.0d0 ! e2(3)=0.0d0 ! endif !fail ! print *,i,'a tu?' ! do j=1,3 ! c(j,i)=c(j,i-1)-1.9d0*e2(j) ! enddo ! else !unres_pdb do j=1,3 dcj=(c(j,i-2)-c(j,i-3))/2.0 if (dcj.eq.0) dcj=1.23591524223 c(j,i)=c(j,i-1)+dcj c(j,nres+i)=c(j,i) enddo ! endif !unres_pdb else !itype(i+1,1).eq.ntyp1 ! if (unres_pdb) then ! 2/15/2013 by Adam: corrected insertion of the first dummy residue ! call refsys(i+1,i+2,i+3,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,i)=c(j,i+1)-1.9d0*e2(j) c(j,i)=c(j,i-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0) enddo ! else !unres_pdb do j=1,3 dcj=(c(j,i+3)-c(j,i+2))/2.0 if (dcj.eq.0) dcj=1.23591524223 c(j,i)=c(j,i+1)-dcj c(j,nres+i)=c(j,i) enddo ! endif !unres_pdb endif !itype(i+1,1).eq.ntyp1 endif !itype.eq.ntyp1 enddo 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 ! print *,"molecule",molecule if ((itype(nres,1).ne.10)) then nres=nres+1 if (molecule.eq.5) molecule=molecprev itype(nres,molecule)=ntyp1_molec(molecule) nres_molec(molecule)=nres_molec(molecule)+1 ! 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)-1.9d0*e2(j) ! enddo ! else do j=1,3 dcj=(c(j,nres-2)-c(j,nres-3))/2.0 c(j,nres)=c(j,nres-1)+dcj c(j,2*nres)=c(j,nres) enddo ! endif endif ! print *,'I have',nres, nres_molec(:) !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,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)-1.9d0*e2(j) c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/sqrt(2.0d0) enddo else do j=1,3 dcj=(c(j,4)-c(j,3))/2.0 c(j,1)=c(j,2)-dcj c(j,nres+1)=c(j,1) enddo endif endif ! First lets assign correct dummy to correct type of chain ! 1) First residue if (itype(1,1).eq.ntyp1) then if (itype(2,1).eq.0) then do j=2,5 if (itype(2,j).ne.0) then itype(1,1)=0 itype(1,j)=ntyp1_molec(j) nres_molec(1)=nres_molec(1)-1 nres_molec(j)=nres_molec(j)+1 go to 3231 endif enddo 3231 continue endif endif print *,'I have',nres, nres_molec(:) ! 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,16x,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,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') & (restyp(itype(ires,j),j),j=1,5),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,'(i5,i3,2x,a,3f8.3,5x,3f8.3)') & ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),& (c(j,nres+ires),j=1,3) enddo endif ! NOW LETS ROCK! SORTING allocate(c_temporary(3,2*nres)) allocate(itype_temporary(nres,5)) if (.not.allocated(molnum)) allocate(molnum(nres+1)) if (.not.allocated(istype)) write(iout,*) & "SOMETHING WRONG WITH ISTYTPE" allocate(istype_temp(nres)) itype_temporary(:,:)=0 seqalingbegin=1 do k=1,5 do i=1,nres if (itype(i,k).ne.0) then do j=1,3 c_temporary(j,seqalingbegin)=c(j,i) c_temporary(j,seqalingbegin+nres)=c(j,i+nres) enddo itype_temporary(seqalingbegin,k)=itype(i,k) print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin istype_temp(seqalingbegin)=istype(i) molnum(seqalingbegin)=k seqalingbegin=seqalingbegin+1 endif enddo enddo do i=1,2*nres do j=1,3 c(j,i)=c_temporary(j,i) enddo enddo do k=1,5 do i=1,nres itype(i,k)=itype_temporary(i,k) istype(i)=istype_temp(i) enddo enddo ! if (itype(1,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)-1.9d0*e2(j) ! enddo ! else ! do j=1,3 ! dcj=(c(j,4)-c(j,3))/2.0 ! c(j,1)=c(j,2)-dcj ! c(j,nres+1)=c(j,1) ! enddo ! endif ! endif if (lprn) then write (iout,'(/a)') & "Cartesian coordinates of the reference structure after sorting" write (iout,'(a,16x,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,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') & (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),& (c(j,ires+nres),j=1,3) enddo endif ! print *,seqalingbegin,nres 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 write (iout,*) "symetr", symetr do i=1,nres lll=lll+1 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3) if (i.gt.1) then if ((itype(i-1,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,1),(chain_rep(j,kkk,i), j=1,3) ! enddo ! enddo ! enddiagnostic ! makes copy of chains write (iout,*) "symetr", symetr do j=1,3 dc(j,0)=c(j,1) enddo 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,1),1),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,'(',i5,')',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!,& integer i 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) !C SHIELD keyword sets if the shielding effect of side-chains is used !C 0 denots no shielding is used all peptide are equally despite the !C solvent accesible area !C 1 the newly introduced function !C 2 reseved for further possible developement call readi(controlcard,'SHIELD',shield_mode,0) !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then write(iout,*) "shield_mode",shield_mode call readi(controlcard,'TORMODE',tor_mode,0) !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then write(iout,*) "torsional and valence angle mode",tor_mode !C Varibles set size of box with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0 protein=index(controlcard,"PROTEIN").gt.0 ions=index(controlcard,"IONS").gt.0 nucleic=index(controlcard,"NUCLEIC").gt.0 write (iout,*) "with_theta_constr ",with_theta_constr AFMlog=(index(controlcard,'AFM')) selfguide=(index(controlcard,'SELFGUIDE')) print *,'AFMlog',AFMlog,selfguide,"KUPA" call readi(controlcard,'GENCONSTR',genconstr,0) call reada(controlcard,'BOXX',boxxsize,100.0d0) call reada(controlcard,'BOXY',boxysize,100.0d0) call reada(controlcard,'BOXZ',boxzsize,100.0d0) call readi(controlcard,'TUBEMOD',tubemode,0) print *,"SCELE",scelemode call readi(controlcard,"SCELEMODE",scelemode,0) print *,"SCELE",scelemode ! elemode = 0 is orignal UNRES electrostatics ! elemode = 1 is "Momo" potentials in progress ! elemode = 2 is in development EVALD write (iout,*) TUBEmode,"TUBEMODE" if (TUBEmode.gt.0) then call reada(controlcard,"XTUBE",tubecenter(1),0.0d0) call reada(controlcard,"YTUBE",tubecenter(2),0.0d0) call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0) call reada(controlcard,"RTUBE",tubeR0,0.0d0) call reada(controlcard,"TUBETOP",bordtubetop,boxzsize) call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0) call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0) buftubebot=bordtubebot+tubebufthick buftubetop=bordtubetop-tubebufthick endif ! CUTOFFF ON ELECTROSTATICS call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0) call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0) write(iout,*) "R_CUT_ELE=",r_cut_ele ! Lipidic parameters call reada(controlcard,"LIPTHICK",lipthick,0.0d0) call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0) if (lipthick.gt.0.0d0) then bordliptop=(boxzsize+lipthick)/2.0 bordlipbot=bordliptop-lipthick if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE" buflipbot=bordlipbot+lipbufthick bufliptop=bordliptop-lipbufthick if ((lipbufthick*2.0d0).gt.lipthick) & write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF" endif !lipthick.gt.0 write(iout,*) "bordliptop=",bordliptop write(iout,*) "bordlipbot=",bordlipbot write(iout,*) "bufliptop=",bufliptop write(iout,*) "buflipbot=",buflipbot write (iout,*) "SHIELD MODE",shield_mode !C------------------------- 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.' if (shield_mode.gt.0) then pi=4.0D0*datan(1.0D0) !C VSolvSphere the volume of solving sphere print *,pi,"pi" !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it !C there will be no distinction between proline peptide group and normal peptide !C group in case of shielding parameters VSolvSphere=4.0/3.0*pi*(4.50d0)**3 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3 write (iout,*) VSolvSphere,VSolvSphere_div !C long axis of side chain ! do i=1,ntyp ! long_r_sidechain(i)=vbldsc0(1,i) ! short_r_sidechain(i)=sigma0(i) ! write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),& ! sigma0(i) ! enddo buff_shield=1.0d0 endif 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,j 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 preminim=(index(controlcard,'PREMINIM').gt.0) write (iout,*) "PREMINIM ",preminim dccart=(index(controlcard,'CART').gt.0) if (preminim) call read_minim ! 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 ! allocate(gamp do j=1,5 !types of molecules gamp(j)=scal_fric*(pstok(j)+rwat)*eta stdfp(j)=dsqrt(2*Rb*t_bath/d_time) enddo allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1) do j=1,5 !types of molecules do i=1,ntyp gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time) enddo 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(1),stdfp*dsqrt(gamp(1)) do i=1,ntyp write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),& gamsc(i,1),stdfsc(i,1)*dsqrt(gamsc(i,1)) 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 MD_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('BONDPAR_NUCL',bondname_nucl) open (ibond_nucl,file=bondname_nucl,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) call getenv_loc('THETPAR_NUCL',thetname_nucl) open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared) call getenv_loc('ROTPAR_NUCL',rotname_nucl) open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared) call getenv_loc('TORPAR_NUCL',torname_nucl) open (itorp_nucl,file=torname_nucl,status='old',readonly,shared) call getenv_loc('TORDPAR_NUCL',tordname_nucl) open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared) call getenv_loc('SIDEPAR_NUCL',sidename_nucl) open (isidep_nucl,file=sidename_nucl,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') call getenv_loc('BONDPAR_NUCL',bondname_nucl) open (ibond_nucl,file=bondname_nucl,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') call getenv_loc('THETPAR_NUCL',thetname_nucl) open (ithep_nucl,file=thetname_nucl,status='old',action='read') call getenv_loc('ROTPAR_NUCL',rotname_nucl) open (irotam_nucl,file=rotname_nucl,status='old',action='read') call getenv_loc('TORPAR_NUCL',torname_nucl) open (itorp_nucl,file=torname_nucl,status='old',action='read') call getenv_loc('TORDPAR_NUCL',tordname_nucl) open (itordp_nucl,file=tordname_nucl,status='old',action='read') call getenv_loc('SIDEPAR_NUCL',sidename_nucl) open (isidep_nucl,file=sidename_nucl,status='old',action='read') call getenv_loc('LIPTRANPAR',liptranname) open (iliptranpar,file=liptranname,status='old',action='read') call getenv_loc('TUBEPAR',tubename) open (itube,file=tubename,status='old',action='read') call getenv_loc('IONPAR',ionname) open (iion,file=ionname,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('BONDPAR_NUCL',bondname_nucl) open (ibond_nucl,file=bondname_nucl,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') open (ithep_nucl,file=thetname_nucl,status='old') call getenv_loc('ROTPAR_NUCL',rotname_nucl) open (irotam_nucl,file=rotname_nucl,status='old') call getenv_loc('TORPAR_NUCL',torname_nucl) open (itorp_nucl,file=torname_nucl,status='old') call getenv_loc('TORDPAR_NUCL',tordname_nucl) open (itordp_nucl,file=tordname_nucl,status='old') call getenv_loc('SIDEPAR_NUCL',sidename_nucl) open (isidep_nucl,file=sidename_nucl,status='old') call getenv_loc('LIPTRANPAR',liptranname) open (iliptranpar,file=liptranname,status='old') call getenv_loc('TUBEPAR',tubename) open (itube,file=tubename,status='old') call getenv_loc('IONPAR',ionname) open (iion,file=ionname,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('BONDPAR_NUCL',bondname_nucl) open (ibond_nucl,file=bondname_nucl,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) call getenv_loc('THETPAR_NUCL',thetname_nucl) open (ithep_nucl,file=thetname_nucl,status='old',action='read') call getenv_loc('ROTPAR_NUCL',rotname_nucl) open (irotam_nucl,file=rotname_nucl,status='old',action='read') call getenv_loc('TORPAR_NUCL',torname_nucl) open (itorp_nucl,file=torname_nucl,status='old',action='read') call getenv_loc('TORDPAR_NUCL',tordname_nucl) open (itordp_nucl,file=tordname_nucl,status='old',action='read') call getenv_loc('SIDEPAR_NUCL',sidename_nucl) open (isidep_nucl,file=sidename_nucl,status='old',action='read') call getenv_loc('SIDEPAR_SCBASE',sidename_scbase) open (isidep_scbase,file=sidename_scbase,status='old',action='read') call getenv_loc('PEPPAR_PEPBASE',pepname_pepbase) open (isidep_pepbase,file=pepname_pepbase,status='old',action='read') call getenv_loc('SCPAR_PHOSPH',pepname_scpho) open (isidep_scpho,file=pepname_scpho,status='old',action='read') call getenv_loc('PEPPAR_PHOSPH',pepname_peppho) open (isidep_peppho,file=pepname_peppho,status='old',action='read') call getenv_loc('LIPTRANPAR',liptranname) open (iliptranpar,file=liptranname,status='old',action='read') call getenv_loc('TUBEPAR',tubename) open (itube,file=tubename,status='old',action='read') call getenv_loc('IONPAR',ionname) open (iion,file=ionname,status='old',action='read') #ifndef CRYST_SC call getenv_loc('ROTPARPDB',rotname_pdb) open (irotam_pdb,file=rotname_pdb,status='old',action='read') #endif #endif call getenv_loc('SCPPAR_NUCL',scpname_nucl) #if defined(WINIFL) || defined(WINPGI) open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared) #elif (defined CRAY) || (defined AIX) open (iscpp_nucl,file=scpname_nucl,status='old',action='read') #elif (defined G77) open (iscpp_nucl,file=scpname_nucl,status='old') #else open (iscpp_nucl,file=scpname_nucl,status='old',action='read') #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 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 totTafm=totT ! do i=1,2*nres ! AL 4/17/17: Now reading d_t(0,:) too do i=0,2*nres read(irest2,'(3e15.5)') (d_t(j,i),j=1,3) enddo ! do i=1,2*nres ! AL 4/17/17: Now reading d_c(0,:) too do i=0,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