2 !-----------------------------------------------------------------------
6 !-----------------------------------------------------------------------------
7 ! Max. number of AA residues
8 integer,parameter :: maxres=6000!1200
9 ! Appr. max. number of interaction sites
10 integer,parameter :: maxres2=2*maxres
11 ! parameter (maxres6=6*maxres)
12 ! parameter (mmaxres2=(maxres2*(maxres2+1)/2))
13 !-----------------------------------------------------------------------------
14 ! Max. number of S-S bridges
15 ! integer,parameter :: maxss=20
16 !-----------------------------------------------------------------------------
17 ! Max. number of derivatives of virtual-bond and side-chain vectors in theta
19 ! integer,parameter :: maxdim=(maxres-1)*(maxres-2)/2
20 !-----------------------------------------------------------------------------
23 !-----------------------------------------------------------------------------
25 !-----------------------------------------------------------------------------
27 !-----------------------------------------------------------------------------
28 subroutine read_bridge
29 ! Read information about disulfide bridges.
30 use geometry_data, only: nres
32 use control_data, only:out1file
34 ! implicit real*8 (a-h,o-z)
35 ! include 'DIMENSIONS'
39 ! include 'COMMON.IOUNITS'
40 ! include 'COMMON.GEO'
41 ! include 'COMMON.VAR'
42 ! include 'COMMON.INTERACT'
43 ! include 'COMMON.LOCAL'
44 ! include 'COMMON.NAMES'
45 ! include 'COMMON.CHAIN'
46 ! include 'COMMON.FFIELD'
47 ! include 'COMMON.SBRIDGE'
48 ! include 'COMMON.HEADER'
49 ! include 'COMMON.CONTROL'
50 ! include 'COMMON.DBASE'
51 ! include 'COMMON.THREAD'
52 ! include 'COMMON.TIME1'
53 ! include 'COMMON.SETUP'
57 ! Read bridging residues.
61 read (inp,*) (iss(i),i=1,ns)
65 if(me.eq.king.or..not.out1file) &
66 write (iout,*) 'ns=',ns
68 write(iout,*) ' iss:',(iss(i),i=1,ns)
69 ! Check whether the specified bridging residues are cystines.
71 if (itype(iss(i)).ne.1) then
72 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)') &
73 'Do you REALLY think that the residue ',&
74 restyp(itype(iss(i))),i,&
75 ' can form a disulfide bridge?!!!'
76 write (*,'(2a,i3,a)') &
77 'Do you REALLY think that the residue ',&
78 restyp(itype(iss(i))),i,&
79 ' can form a disulfide bridge?!!!'
81 call MPI_Finalize(MPI_COMM_WORLD,ierror)
86 ! Read preformed bridges.
90 if(.not.allocated(ihpb)) allocate(ihpb(nss))
91 if(.not.allocated(jhpb)) allocate(jhpb(nss))
92 read (inp,*) (ihpb(i),jhpb(i),i=1,nss)
94 if(.not.allocated(dhpb)) allocate(dhpb(nss))
95 if(.not.allocated(forcon)) allocate(forcon(nss))!(maxdim) !el maxdim=(maxres-1)*(maxres-2)/2
96 if(.not.allocated(dhpb1)) allocate(dhpb1(nss))
97 if(.not.allocated(ibecarb)) allocate(ibecarb(nss))
98 ! el Initialize the bridge array
103 !--------------------
105 write(iout,*)'nss=',nss !el,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
107 write(iout,*)'ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
109 ! Check if the residues involved in bridges are in the specified list of
113 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j) &
114 .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
115 write (iout,'(a,i3,a)') 'Disulfide pair',i,&
116 ' contains residues present in other pairs.'
117 write (*,'(a,i3,a)') 'Disulfide pair',i,&
118 ' contains residues present in other pairs.'
120 call MPI_Finalize(MPI_COMM_WORLD,ierror)
126 if (ihpb(i).eq.iss(j)) goto 10
128 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
131 if (jhpb(i).eq.iss(j)) goto 20
133 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
144 ! write(iout,*) "end read_bridge"
146 end subroutine read_bridge
147 !-----------------------------------------------------------------------------
148 subroutine read_x(kanal,*)
152 use geometry, only:int_from_cart1
153 ! implicit real*8 (a-h,o-z)
154 ! include 'DIMENSIONS'
155 ! include 'COMMON.GEO'
156 ! include 'COMMON.VAR'
157 ! include 'COMMON.CHAIN'
158 ! include 'COMMON.IOUNITS'
159 ! include 'COMMON.CONTROL'
160 ! include 'COMMON.LOCAL'
161 ! include 'COMMON.INTERACT'
162 ! Read coordinates from input
165 integer :: l,k,j,i,kanal
167 read(kanal,'(8f10.5)',end=10,err=10) &
168 ((c(l,k),l=1,3),k=1,nres),&
169 ((c(l,k+nres),l=1,3),k=nnt,nct)
172 c(j,2*nres)=c(j,nres)
174 call int_from_cart1(.false.)
177 dc(j,i)=c(j,i+1)-c(j,i)
178 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
182 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
184 dc(j,i+nres)=c(j,i+nres)-c(j,i)
185 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
192 end subroutine read_x
193 !-----------------------------------------------------------------------------
194 subroutine read_threadbase
199 ! implicit real*8 (a-h,o-z)
200 ! include 'DIMENSIONS'
201 ! include 'COMMON.IOUNITS'
202 ! include 'COMMON.GEO'
203 ! include 'COMMON.VAR'
204 ! include 'COMMON.INTERACT'
205 ! include 'COMMON.LOCAL'
206 ! include 'COMMON.NAMES'
207 ! include 'COMMON.CHAIN'
208 ! include 'COMMON.FFIELD'
209 ! include 'COMMON.SBRIDGE'
210 ! include 'COMMON.HEADER'
211 ! include 'COMMON.CONTROL'
212 ! include 'COMMON.DBASE'
213 ! include 'COMMON.THREAD'
214 ! include 'COMMON.TIME1'
219 ! Read pattern database for threading.
221 allocate(cart_base(3,maxres_base,nseq)) !(3,maxres_base,maxseq)
222 allocate(nres_base(3,nseq)) !(3,maxseq)
223 allocate(str_nam(nseq)) !(maxseq)
225 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),&
226 nres_base(2,i),nres_base(3,i)
227 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,&
229 ! write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
230 ! & nres_base(2,i),nres_base(3,i)
231 ! write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
235 if (weidis.eq.0.0D0) weidis=0.1D0
244 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
245 write (iout,'(a,i5)') 'nexcl: ',nexcl
246 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
248 end subroutine read_threadbase
249 !-----------------------------------------------------------------------------
251 !el subroutine read_angles(kanal,iscor,energ,iprot,*)
252 subroutine read_angles(kanal,*)
256 ! implicit real*8 (a-h,o-z)
257 ! include 'DIMENSIONS'
258 ! include 'DIMENSIONS.ZSCOPT'
259 ! include 'COMMON.INTERACT'
260 ! include 'COMMON.SBRIDGE'
261 ! include 'COMMON.GEO'
262 ! include 'COMMON.VAR'
263 ! include 'COMMON.CHAIN'
264 ! include 'COMMON.IOUNITS'
265 character(len=80) :: lineh
266 integer :: iscor,iprot,ic
267 real(kind=8) :: energ
269 subroutine read_angles(kanal,*)
274 ! implicit real*8 (a-h,o-z)
275 ! include 'DIMENSIONS'
276 ! include 'COMMON.GEO'
277 ! include 'COMMON.VAR'
278 ! include 'COMMON.CHAIN'
279 ! include 'COMMON.IOUNITS'
280 ! include 'COMMON.CONTROL'
282 ! Read angles from input
287 read(kanal,'(a80)',end=10,err=10) lineh
288 read(lineh(:5),*,err=8) ic
289 read(lineh(6:),*,err=8) energ
292 print *,'error, assuming e=1d10',lineh
296 read(lineh(18:),*,end=10,err=10) nss
298 read (lineh(20:),*,end=10,err=10) &
299 (IHPB(I),JHPB(I),I=1,NSS),iscor
301 read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
302 read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),&
305 ! print *,"energy",energ," iscor",iscor
307 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
308 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
309 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
310 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
313 ! 9/7/01 avoid 180 deg valence angle
314 if (theta(i).gt.179.99d0) theta(i)=179.99d0
316 theta(i)=deg2rad*theta(i)
317 phi(i)=deg2rad*phi(i)
318 alph(i)=deg2rad*alph(i)
319 omeg(i)=deg2rad*omeg(i)
323 end subroutine read_angles
324 !-----------------------------------------------------------------------------
325 subroutine reada(rekord,lancuch,wartosc,default)
328 character*(*) :: rekord,lancuch
329 real(kind=8) :: wartosc,default
330 integer :: iread !,ilen
332 iread=index(rekord,lancuch)
337 iread=iread+ilen(lancuch)+1
338 read (rekord(iread:),*,err=10,end=10) wartosc
343 !-----------------------------------------------------------------------------
344 subroutine readi(rekord,lancuch,wartosc,default)
347 character*(*) :: rekord,lancuch
348 integer :: wartosc,default
349 integer :: iread !,ilen
351 iread=index(rekord,lancuch)
356 iread=iread+ilen(lancuch)+1
357 read (rekord(iread:),*,err=10,end=10) wartosc
362 !-----------------------------------------------------------------------------
363 subroutine multreadi(rekord,lancuch,tablica,dim,default)
367 integer :: tablica(dim),default
368 character*(*) :: rekord,lancuch
369 character(len=80) :: aux
370 integer :: iread !,ilen
375 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
376 if (iread.eq.0) return
377 iread=iread+ilen(lancuch)+1
378 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
380 end subroutine multreadi
381 !-----------------------------------------------------------------------------
382 subroutine multreada(rekord,lancuch,tablica,dim,default)
386 real(kind=8) :: tablica(dim),default
387 character*(*) :: rekord,lancuch
388 character(len=80) :: aux
389 integer :: iread !,ilen
394 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
395 if (iread.eq.0) return
396 iread=iread+ilen(lancuch)+1
397 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
399 end subroutine multreada
400 !-----------------------------------------------------------------------------
401 subroutine card_concat(card,to_upper)
403 ! dla UNRESA to_upper jest zawsze .true.
404 ! implicit real*8 (a-h,o-z)
405 ! include 'DIMENSIONS'
406 ! include 'COMMON.IOUNITS'
408 character(len=80) :: karta !,ucase
411 read (inp,'(a)') karta
412 if (to_upper) karta=ucase(karta)
414 do while (karta(80:80).eq.'&')
415 card=card(:ilen(card)+1)//karta(:79)
416 read (inp,'(a)') karta
417 if (to_upper) karta=ucase(karta)
419 card=card(:ilen(card)+1)//karta
421 end subroutine card_concat
422 !-----------------------------------------------------------------------------
423 subroutine read_dist_constr
426 use geometry, only: dist
430 ! implicit real*8 (a-h,o-z)
431 ! include 'DIMENSIONS'
435 ! include 'COMMON.SETUP'
436 ! include 'COMMON.CONTROL'
437 ! include 'COMMON.CHAIN'
438 ! include 'COMMON.IOUNITS'
439 ! include 'COMMON.SBRIDGE'
440 integer,dimension(2,100) :: ifrag_,ipair_
441 real(kind=8),dimension(100) :: wfrag_,wpair_
442 character(len=640) :: controlcard
445 integer :: i,k,j,ddjk,ii,jj,itemp
446 integer :: nfrag_,npair_,ndist_
447 real(kind=8) :: dist_cut
449 ! write (iout,*) "Calling read_dist_constr"
450 ! write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
452 call card_concat(controlcard,.true.)
453 call readi(controlcard,"NFRAG",nfrag_,0)
454 call readi(controlcard,"NPAIR",npair_,0)
455 call readi(controlcard,"NDIST",ndist_,0)
456 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
457 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
458 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
459 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
460 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
461 ! write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
462 ! write (iout,*) "IFRAG"
464 ! write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
466 ! write (iout,*) "IPAIR"
468 ! write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
470 if(.not.allocated(ihpb)) allocate(ihpb(maxdim))
471 if(.not.allocated(jhpb)) allocate(jhpb(maxdim))
472 if(.not.allocated(dhpb)) allocate(dhpb(maxdim))
473 if(.not.allocated(forcon)) allocate(forcon(maxdim))
477 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
478 if (ifrag_(2,i).gt.nstart_sup+nsup-1) &
479 ifrag_(2,i)=nstart_sup+nsup-1
480 ! write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
482 if (wfrag_(i).gt.0.0d0) then
483 do j=ifrag_(1,i),ifrag_(2,i)-1
485 ! write (iout,*) "j",j," k",k
487 if (constr_dist.eq.1) then
492 forcon(nhpb)=wfrag_(i)
493 else if (constr_dist.eq.2) then
494 if (ddjk.le.dist_cut) then
499 forcon(nhpb)=wfrag_(i)
506 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
509 if (.not.out1file .or. me.eq.king) &
510 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",&
511 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
513 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",&
514 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
521 if (wpair_(i).gt.0.0d0) then
529 do j=ifrag_(1,ii),ifrag_(2,ii)
530 do k=ifrag_(1,jj),ifrag_(2,jj)
534 forcon(nhpb)=wpair_(i)
537 if (.not.out1file .or. me.eq.king) &
538 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
539 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
541 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
542 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
549 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
550 if (forcon(nhpb+1).gt.0.0d0) then
552 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
554 if (.not.out1file .or. me.eq.king) &
555 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
556 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
558 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
559 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
565 end subroutine read_dist_constr
566 !-----------------------------------------------------------------------------
578 !-----------------------------------------------------------------------------
579 subroutine copy_to_tmp(source)
581 ! include "DIMENSIONS"
582 ! include "COMMON.IOUNITS"
583 character*(*) :: source
584 character(len=256) :: tmpfile
588 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
589 inquire(file=tmpfile,exist=ex)
591 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),&
592 " to temporary directory..."
593 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
594 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
597 end subroutine copy_to_tmp
598 !-----------------------------------------------------------------------------
599 subroutine move_from_tmp(source)
601 ! include "DIMENSIONS"
602 ! include "COMMON.IOUNITS"
603 character*(*) :: source
606 write (*,*) "Moving ",source(:ilen(source)),&
607 " from temporary directory to working directory"
608 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
609 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
611 end subroutine move_from_tmp
612 !-----------------------------------------------------------------------------
614 !-----------------------------------------------------------------------------
615 ! $Date: 1994/10/12 17:24:21 $
618 logical function find_arg(ipos,line,errflag)
620 integer, parameter :: maxlen = 80
621 character(len=80) :: line
622 character(len=1) :: empty=' ',equal='='
625 ! This function returns .TRUE., if an argument follows keyword keywd; if so
626 ! IPOS will point to the first non-blank character of the argument. Returns
627 ! .FALSE., if no argument follows the keyword; in this case IPOS points
628 ! to the first non-blank character of the next keyword.
630 do while (line(ipos:ipos) .eq. empty .and. ipos.le.maxlen)
634 if (line(ipos:ipos).eq.equal) then
637 do while (line(ipos:ipos) .eq. empty .and. ipos.le.maxlen)
640 if (ipos.gt.maxlen) errflag=.true.
646 end function find_arg
647 !-----------------------------------------------------------------------------
648 logical function find_group(iunit,jout,key1)
650 character*(*) :: key1
651 character(len=80) :: karta !,ucase
652 integer :: iunit,jout
659 do while (index(ucase(karta),key1(1:ll)).eq.0.or.lcom(1,karta))
660 read (iunit,'(a)',end=10) karta
662 write (jout,'(2a)') '> ',karta(1:78)
665 10 find_group=.false.
667 end function find_group
668 !-----------------------------------------------------------------------------
669 logical function iblnk(charc)
670 character(len=1) :: charc
673 iblnk = (n.eq.9) .or. (n.eq.10) .or. (charc.eq.' ')
676 !-----------------------------------------------------------------------------
677 integer function ilen(string)
678 character*(*) :: string
682 1 if ( ilen .gt. 0 ) then
683 if ( iblnk( string(ilen:ilen) ) ) then
690 !-----------------------------------------------------------------------------
691 integer function in_keywd_set(nkey,ikey,narg,keywd,keywdset)
692 integer :: nkey,i,ikey,narg
693 character(len=16) :: keywd,keywdset(1:nkey,0:nkey)
694 ! character(len=16) :: ucase
697 if (ucase(keywd).eq.keywdset(i,ikey)) then
703 ! No match to the allowed set of keywords if this point is reached.
706 end function in_keywd_set
707 !-----------------------------------------------------------------------------
708 character function lcase(string)
709 integer :: i, k, idiff
710 character*(*) :: string
711 character(len=1) :: c
712 character(len=40) :: chtmp
718 if (string(k+1:) .ne. ' ') then
722 idiff = ichar('a') - ichar('A')
726 if (lge(c,'A') .and. lle(c,'Z')) then
727 lcase(i:i) = char(ichar(c) + idiff)
732 !-----------------------------------------------------------------------------
733 logical function lcom(ipos,karta)
734 character(len=80) :: karta
735 character :: koment(2) = (/'!','#'/)
740 if (karta(ipos:ipos).eq.koment(i)) lcom=.true.
744 !-----------------------------------------------------------------------------
745 logical function lower_case(ch)
747 lower_case=(ch.ge.'a' .and. ch.le.'z')
749 end function lower_case
750 !-----------------------------------------------------------------------------
751 subroutine mykey(line,keywd,ipos,blankline,errflag)
753 ! This subroutine seeks a non-empty substring keywd in the string LINE.
754 ! The substring begins with the first character different from blank and
755 ! "=" encountered right to the pointer IPOS (inclusively) and terminates
756 ! at the character left to the first blank or "=". When the subroutine is
757 ! exited, the pointer IPOS is moved to the position of the terminator in LINE.
758 ! The logical variable BLANKLINE is set at .TRUE., if LINE(IPOS:) contains
759 ! only separators or the maximum length of the data line (80) has been reached.
760 ! The logical variable ERRFLAG is set at .TRUE. if the string
761 ! consists only from a "=".
762 integer, parameter :: maxlen=80
763 character(len=1) :: empty=' ',equal='=',comma=','
764 character*(*) :: keywd
765 character(len=80) :: line
766 logical :: blankline,errflag !EL,lcom
767 integer :: ipos,istart,iend
769 do while (line(ipos:ipos).eq.empty .and. (ipos.le.maxlen))
772 if (ipos.gt.maxlen .or. lcom(ipos,line) ) then
773 ! At this point the rest of the input line turned out to contain only blanks
774 ! or to be commented out.
780 ! Checks whether the current char is a separator.
781 do while (line(ipos:ipos).ne.empty .and. line(ipos:ipos).ne.equal &
782 .and. line(ipos:ipos).ne.comma .and. ipos.le.maxlen)
786 ! Error flag set to .true., if the length of the keyword was found less than 1.
787 if (iend.lt.istart) then
791 keywd=line(istart:iend)
794 !-----------------------------------------------------------------------------
795 subroutine numstr(inum,numm)
796 character(len=10) :: huj='0123456789'
797 character*(*) :: numm
798 integer :: inumm,inum,inum1,inum2
803 numm(3:3)=huj(inum2+1:inum2+1)
807 numm(2:2)=huj(inum2+1:inum2+1)
811 numm(1:1)=huj(inum2+1:inum2+1)
813 end subroutine numstr
814 !-----------------------------------------------------------------------------
815 function ucase(string)
816 integer :: i, k, idiff
817 character(*) :: string
818 character(len=len(string)) :: ucase
819 character(len=1) :: c
820 character(len=40) :: chtmp
826 if (string(k+1:) .ne. ' ') then
830 idiff = ichar('a') - ichar('A')
834 if (lge(c,'a') .and. lle(c,'z')) then
835 ucase(i:i) = char(ichar(c) - idiff)
840 !-----------------------------------------------------------------------------
842 !-----------------------------------------------------------------------------
843 subroutine pdbout(etot,tytul,iunit)
845 use geometry_data, only: c,nres
850 ! implicit real*8 (a-h,o-z)
851 ! include 'DIMENSIONS'
852 ! include 'COMMON.CHAIN'
853 ! include 'COMMON.INTERACT'
854 ! include 'COMMON.NAMES'
855 ! include 'COMMON.IOUNITS'
856 ! include 'COMMON.HEADER'
857 ! include 'COMMON.SBRIDGE'
858 ! include 'COMMON.DISTFIT'
859 ! include 'COMMON.MD'
860 !el character(len=50) :: tytul
861 character*(*) :: tytul
862 character(len=1),dimension(10) :: chainid= (/'A','B','C','D','E','F','G','H','I','J'/)
863 integer,dimension(nres) :: ica !(maxres)
866 integer :: j,iti,itj,itk,itl,i,iatom,ichain,ires,iunit
871 if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
873 write (iunit,'(3a,1pe15.5)') 'REMARK ',tytul,' ENERGY ',etot
874 !model write (iunit,'(a5,i6)') 'MODEL',1
875 if (nhfrag.gt.0) then
877 iti=itype(hfrag(1,j))
878 itj=itype(hfrag(2,j))
880 write (iunit,'(a5,i5,1x,a1,i1,2x,a3,i7,2x,a3,i7,i3,t76,i5)') &
882 restyp(iti),hfrag(1,j)-1,&
883 restyp(itj),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j)
885 write (iunit,'(a5,i5,1x,a1,i2,1x,a3,i7,2x,a3,i7,i3)') &
887 restyp(iti),hfrag(1,j)-1,&
888 restyp(itj),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j)
893 if (nbfrag.gt.0) then
897 iti=itype(bfrag(1,j))
898 itj=itype(bfrag(2,j)-1)
900 write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3)') &
902 restyp(iti),bfrag(1,j)-1,&
903 restyp(itj),bfrag(2,j)-2,0
905 if (bfrag(3,j).gt.bfrag(4,j)) then
907 itk=itype(bfrag(3,j))
908 itl=itype(bfrag(4,j)+1)
910 write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3,2x,a1,2x,a3,i6,3x,a1,2x,a3,i6)') &
912 restyp(itl),bfrag(4,j),&
913 restyp(itk),bfrag(3,j)-1,-1,&
914 "N",restyp(itk),bfrag(3,j)-1,&
915 "O",restyp(iti),bfrag(1,j)-1
919 itk=itype(bfrag(3,j))
920 itl=itype(bfrag(4,j)-1)
923 write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3,2x,a1,2x,a3,i6,3x,a1,2x,a3,i6)') &
925 restyp(itk),bfrag(3,j)-1,&
926 restyp(itl),bfrag(4,j)-2,1,&
927 "N",restyp(itk),bfrag(3,j)-1,&
928 "O",restyp(iti),bfrag(1,j)-1
940 write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') &
941 'SSBOND',i,'CYS',idssb(i)-nnt+1,&
944 write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') &
945 'SSBOND',i,'CYS',ihpb(i)-nnt+1-nres,&
946 'CYS',jhpb(i)-nnt+1-nres
957 if ((iti.eq.ntyp1).and.(iti1.eq.ntyp1)) cycle
958 if (iti.eq.ntyp1) then
961 write (iunit,'(a)') 'TER'
966 write (iunit,10) iatom,restyp(iti),chainid(ichain),&
967 ires,(c(j,i),j=1,3),vtot(i)
970 write (iunit,20) iatom,restyp(iti),chainid(ichain),&
971 ires,(c(j,nres+i),j=1,3),&
976 write (iunit,'(a)') 'TER'
978 if (itype(i).eq.ntyp1) cycle
979 if (itype(i).eq.10 .and. itype(i+1).ne.ntyp1) then
980 write (iunit,30) ica(i),ica(i+1)
981 else if (itype(i).ne.10 .and. itype(i+1).ne.ntyp1) then
982 write (iunit,30) ica(i),ica(i+1),ica(i)+1
983 else if (itype(i).ne.10 .and. itype(i+1).eq.ntyp1) then
984 write (iunit,30) ica(i),ica(i)+1
987 if (itype(nct).ne.10) then
988 write (iunit,30) ica(nct),ica(nct)+1
992 write (iunit,30) ica(idssb(i))+1,ica(jdssb(i))+1
994 write (iunit,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
997 write (iunit,'(a6)') 'ENDMDL'
998 10 FORMAT ('ATOM',I7,' CA ',A3,1X,A1,I4,4X,3F8.3,f15.3)
999 20 FORMAT ('ATOM',I7,' CB ',A3,1X,A1,I4,4X,3F8.3,f15.3)
1000 30 FORMAT ('CONECT',8I5)
1002 end subroutine pdbout
1003 !-----------------------------------------------------------------------------
1004 subroutine MOL2out(etot,tytul)
1005 ! Prints the Cartesian coordinates of the alpha-carbons in the Tripos mol2
1007 use geometry_data, only: c
1010 ! implicit real*8 (a-h,o-z)
1011 ! include 'DIMENSIONS'
1012 ! include 'COMMON.CHAIN'
1013 ! include 'COMMON.INTERACT'
1014 ! include 'COMMON.NAMES'
1015 ! include 'COMMON.IOUNITS'
1016 ! include 'COMMON.HEADER'
1017 ! include 'COMMON.SBRIDGE'
1018 character(len=32) :: tytul,fd
1019 character(len=3) :: zahl
1020 character(len=6) :: res_num,pom !,ucase
1024 real(kind=8) :: etot
1028 #elif (defined CRAY)
1033 write (imol2,'(a)') '#'
1034 write (imol2,'(a)') &
1035 '# Creating user name: unres'
1036 write (imol2,'(2a)') '# Creation time: ',&
1038 write (imol2,'(/a)') '\@<TRIPOS>MOLECULE'
1039 write (imol2,'(a)') tytul
1040 write (imol2,'(5i5)') nct-nnt+1,nct-nnt+nss+1,nct-nnt+nss+1,0,0
1041 write (imol2,'(a)') 'SMALL'
1042 write (imol2,'(a)') 'USER_CHARGES'
1043 write (imol2,'(a)') '\@<TRIPOS>ATOM'
1045 write (zahl,'(i3)') i
1046 pom=ucase(restyp(itype(i)))
1047 res_num = pom(:3)//zahl(2:)
1048 write (imol2,10) i,(c(j,i),j=1,3),i,res_num,0.0
1050 write (imol2,'(a)') '\@<TRIPOS>BOND'
1052 write (imol2,'(i5,2i6,i2)') i-nnt+1,i-nnt+1,i-nnt+2,1
1055 write (imol2,'(i5,2i6,i2)') nct-nnt+i,ihpb(i),jhpb(i),1
1057 write (imol2,'(a)') '\@<TRIPOS>SUBSTRUCTURE'
1059 write (zahl,'(i3)') i
1060 pom = ucase(restyp(itype(i)))
1061 res_num = pom(:3)//zahl(2:)
1062 write (imol2,30) i-nnt+1,res_num,i-nnt+1,0
1064 10 FORMAT (I7,' CA ',3F10.4,' C.3',I8,1X,A,F11.4,' ****')
1065 30 FORMAT (I7,1x,A,I14,' RESIDUE',I13,' **** ****')
1067 end subroutine MOL2out
1068 !-----------------------------------------------------------------------------
1072 use energy_data, only: itype
1074 ! implicit real*8 (a-h,o-z)
1075 ! include 'DIMENSIONS'
1076 ! include 'COMMON.IOUNITS'
1077 ! include 'COMMON.CHAIN'
1078 ! include 'COMMON.VAR'
1079 ! include 'COMMON.LOCAL'
1080 ! include 'COMMON.INTERACT'
1081 ! include 'COMMON.NAMES'
1082 ! include 'COMMON.GEO'
1083 ! include 'COMMON.TORSION'
1087 write (iout,'(/a)') 'Geometry of the virtual chain.'
1088 write (iout,'(7a)') ' Res ',' d',' Theta',&
1089 ' Phi',' Dsc',' Alpha',' Omega'
1092 write (iout,'(a3,i4,6f10.3)') restyp(iti),i,vbld(i),&
1093 rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),rad2deg*alph(i),&
1097 end subroutine intout
1098 !-----------------------------------------------------------------------------
1100 subroutine briefout(it,ener,free)!,plik)
1102 subroutine briefout(it,ener)
1107 ! implicit real*8 (a-h,o-z)
1108 ! include 'DIMENSIONS'
1109 ! include 'COMMON.IOUNITS'
1110 ! include 'COMMON.CHAIN'
1111 ! include 'COMMON.VAR'
1112 ! include 'COMMON.LOCAL'
1113 ! include 'COMMON.INTERACT'
1114 ! include 'COMMON.NAMES'
1115 ! include 'COMMON.GEO'
1116 ! include 'COMMON.SBRIDGE'
1117 ! print '(a,i5)',intname,igeom
1120 real(kind=8) :: ener,free
1121 ! character(len=80) :: plik
1124 #if defined(AIX) || defined(PGI)
1125 open (igeom,file=intname,position='append')
1127 open (igeom,file=intname,access='append')
1135 WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,NSS)
1137 WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,8)
1139 WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,NSS)
1141 WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,9)
1143 WRITE (igeom,190) (IHPB(I),JHPB(I),I=10,NSS)
1145 ! IF (nvar.gt.nphi) WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
1146 WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
1147 WRITE (igeom,200) (RAD2DEG*PHI(I),I=4,NRES)
1148 ! if (nvar.gt.nphi+ntheta) then
1149 write (igeom,200) (rad2deg*alph(i),i=2,nres-1)
1150 write (igeom,200) (rad2deg*omeg(i),i=2,nres-1)
1153 180 format (I5,F12.3,I2,9(1X,2I3))
1154 190 format (3X,11(1X,2I3))
1157 end subroutine briefout
1158 !-----------------------------------------------------------------------------
1160 subroutine fdate(fd)
1161 character(len=32) :: fd
1164 end subroutine fdate
1166 !-----------------------------------------------------------------------------
1168 real(kind=8) function gyrate(jcon)
1170 real(kind=8) function gyrate()
1173 use geometry_data, only: c
1176 ! implicit real*8 (a-h,o-z)
1177 ! include 'DIMENSIONS'
1178 ! include 'COMMON.INTERACT'
1179 ! include 'COMMON.CHAIN'
1181 real(kind=8),dimension(3) :: cen
1191 cen(j)=cen(j)+c(j,i)
1195 cen(j)=cen(j)/dble(nct-nnt+1)
1200 rg = rg + (c(j,i)-cen(j))**2
1204 gyrate = dsqrt(rg/dble(nct-nnt+1))
1206 gyrate = sqrt(rg/dble(nct-nnt+1))
1211 !-----------------------------------------------------------------------------
1213 subroutine reads(rekord,lancuch,wartosc,default)
1215 character*(*) :: rekord,lancuch,wartosc,default
1216 character(len=80) :: aux
1217 integer :: lenlan,lenrec,iread,ireade
1221 lenlan=ilen(lancuch)
1223 iread=index(rekord,lancuch(:lenlan)//"=")
1224 ! print *,"rekord",rekord," lancuch",lancuch
1225 ! print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
1226 if (iread.eq.0) then
1230 iread=iread+lenlan+1
1231 ! print *,"iread",iread
1232 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
1233 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
1235 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
1237 ! print *,"iread",iread
1238 if (iread.gt.lenrec) then
1243 ! print *,"ireade",ireade
1244 do while (ireade.lt.lenrec .and. &
1245 .not.iblnk(rekord(ireade:ireade)))
1248 wartosc=rekord(iread:ireade)
1250 end subroutine reads
1252 !-----------------------------------------------------------------------------
1254 !-----------------------------------------------------------------------------
1255 subroutine permut(isym)
1257 use geometry_data, only: tabperm
1258 ! implicit real*8 (a-h,o-z)
1259 ! include 'DIMENSIONS'
1260 ! include 'COMMON.LOCAL'
1261 ! include 'COMMON.VAR'
1262 ! include 'COMMON.CHAIN'
1263 ! include 'COMMON.INTERACT'
1264 ! include 'COMMON.IOUNITS'
1265 ! include 'COMMON.GEO'
1266 ! include 'COMMON.NAMES'
1267 ! include 'COMMON.CONTROL'
1272 integer,dimension(isym) :: a
1273 ! parameter(n=symetr)
1286 10 print *,(a(i),i=1,n)
1290 ! write (iout,*) "tututu", kkk
1292 if(nextp(n,a)) go to 10
1294 end subroutine permut
1295 !-----------------------------------------------------------------------------
1296 logical function nextp(n,a)
1298 integer :: n,i,j,k,t
1300 integer,dimension(n) :: a
1302 10 if(a(i).lt.a(i+1)) go to 20
1319 if(a(j).lt.a(i)) go to 40
1326 !-----------------------------------------------------------------------------
1327 !-----------------------------------------------------------------------------