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.
62 read (inp,*) (iss(i),i=1,ns)
66 if(me.eq.king.or..not.out1file) &
67 write (iout,*) 'ns=',ns
69 write(iout,*) ' iss:',(iss(i),i=1,ns)
70 ! Check whether the specified bridging residues are cystines.
72 write(iout,*) i,iss(i)
73 if (itype(iss(i),1).ne.1) then
74 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)') &
75 'Do you REALLY think that the residue ',&
76 restyp(itype(iss(i),1),1),i,&
77 ' can form a disulfide bridge?!!!'
78 write (*,'(2a,i3,a)') &
79 'Do you REALLY think that the residue ',&
80 restyp(itype(iss(i),1),1),i,&
81 ' can form a disulfide bridge?!!!'
83 call MPI_Finalize(MPI_COMM_WORLD,ierror)
89 if(.not.allocated(ihpb)) allocate(ihpb(ns))
90 if(.not.allocated(jhpb)) allocate(jhpb(ns))
92 ! Read preformed bridges.
96 if(.not.allocated(ihpb)) allocate(ihpb(nss))
97 if(.not.allocated(jhpb)) allocate(jhpb(nss))
98 read (inp,*) (ihpb(i),jhpb(i),i=1,nss)
100 if(.not.allocated(dhpb)) allocate(dhpb(nss))
101 if(.not.allocated(forcon)) allocate(forcon(nss))!(maxdim) !el maxdim=(maxres-1)*(maxres-2)/2
102 if(.not.allocated(dhpb1)) allocate(dhpb1(nss))
103 if(.not.allocated(ibecarb)) allocate(ibecarb(nss))
104 ! el Initialize the bridge array
109 !--------------------
111 write(iout,*)'nss=',nss !el,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
113 write(iout,*)'ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
115 ! Check if the residues involved in bridges are in the specified list of
119 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j) &
120 .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
121 write (iout,'(a,i3,a)') 'Disulfide pair',i,&
122 ' contains residues present in other pairs.'
123 write (*,'(a,i3,a)') 'Disulfide pair',i,&
124 ' contains residues present in other pairs.'
126 call MPI_Finalize(MPI_COMM_WORLD,ierror)
132 if (ihpb(i).eq.iss(j)) goto 10
134 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
137 if (jhpb(i).eq.iss(j)) goto 20
139 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
150 write(iout,*) "end read_bridge"
152 end subroutine read_bridge
153 !-----------------------------------------------------------------------------
154 subroutine read_x(kanal,*)
158 use geometry, only:int_from_cart1
159 ! implicit real*8 (a-h,o-z)
160 ! include 'DIMENSIONS'
161 ! include 'COMMON.GEO'
162 ! include 'COMMON.VAR'
163 ! include 'COMMON.CHAIN'
164 ! include 'COMMON.IOUNITS'
165 ! include 'COMMON.CONTROL'
166 ! include 'COMMON.LOCAL'
167 ! include 'COMMON.INTERACT'
168 ! Read coordinates from input
171 integer :: l,k,j,i,kanal
173 read(kanal,'(8f10.5)',end=10,err=10) &
174 ((c(l,k),l=1,3),k=1,nres),&
175 ((c(l,k+nres),l=1,3),k=nnt,nct)
178 c(j,2*nres)=c(j,nres)
180 call int_from_cart1(.false.)
183 dc(j,i)=c(j,i+1)-c(j,i)
184 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
188 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
190 dc(j,i+nres)=c(j,i+nres)-c(j,i)
191 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
198 end subroutine read_x
199 !-----------------------------------------------------------------------------
200 subroutine read_threadbase
205 ! implicit real*8 (a-h,o-z)
206 ! include 'DIMENSIONS'
207 ! include 'COMMON.IOUNITS'
208 ! include 'COMMON.GEO'
209 ! include 'COMMON.VAR'
210 ! include 'COMMON.INTERACT'
211 ! include 'COMMON.LOCAL'
212 ! include 'COMMON.NAMES'
213 ! include 'COMMON.CHAIN'
214 ! include 'COMMON.FFIELD'
215 ! include 'COMMON.SBRIDGE'
216 ! include 'COMMON.HEADER'
217 ! include 'COMMON.CONTROL'
218 ! include 'COMMON.DBASE'
219 ! include 'COMMON.THREAD'
220 ! include 'COMMON.TIME1'
225 ! Read pattern database for threading.
227 allocate(cart_base(3,maxres_base,nseq)) !(3,maxres_base,maxseq)
228 allocate(nres_base(3,nseq)) !(3,maxseq)
229 allocate(str_nam(nseq)) !(maxseq)
231 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),&
232 nres_base(2,i),nres_base(3,i)
233 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,&
235 ! write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
236 ! & nres_base(2,i),nres_base(3,i)
237 ! write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
241 if (weidis.eq.0.0D0) weidis=0.1D0
250 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
251 write (iout,'(a,i5)') 'nexcl: ',nexcl
252 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
254 end subroutine read_threadbase
255 !-----------------------------------------------------------------------------
257 !el subroutine read_angles(kanal,iscor,energ,iprot,*)
258 subroutine read_angles(kanal,*)
262 ! implicit real*8 (a-h,o-z)
263 ! include 'DIMENSIONS'
264 ! include 'DIMENSIONS.ZSCOPT'
265 ! include 'COMMON.INTERACT'
266 ! include 'COMMON.SBRIDGE'
267 ! include 'COMMON.GEO'
268 ! include 'COMMON.VAR'
269 ! include 'COMMON.CHAIN'
270 ! include 'COMMON.IOUNITS'
271 character(len=80) :: lineh
272 integer :: iscor,iprot,ic
273 real(kind=8) :: energ
275 subroutine read_angles(kanal,*)
280 ! implicit real*8 (a-h,o-z)
281 ! include 'DIMENSIONS'
282 ! include 'COMMON.GEO'
283 ! include 'COMMON.VAR'
284 ! include 'COMMON.CHAIN'
285 ! include 'COMMON.IOUNITS'
286 ! include 'COMMON.CONTROL'
288 ! Read angles from input
293 read(kanal,'(a80)',end=10,err=10) lineh
294 read(lineh(:5),*,err=8) ic
295 read(lineh(6:),*,err=8) energ
298 print *,'error, assuming e=1d10',lineh
302 read(lineh(18:),*,end=10,err=10) nss
304 read (lineh(20:),*,end=10,err=10) &
305 (IHPB(I),JHPB(I),I=1,NSS),iscor
307 read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
308 read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),&
311 ! print *,"energy",energ," iscor",iscor
313 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
314 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
315 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
316 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
319 ! 9/7/01 avoid 180 deg valence angle
320 if (theta(i).gt.179.99d0) theta(i)=179.99d0
322 theta(i)=deg2rad*theta(i)
323 phi(i)=deg2rad*phi(i)
324 alph(i)=deg2rad*alph(i)
325 omeg(i)=deg2rad*omeg(i)
329 end subroutine read_angles
330 !-----------------------------------------------------------------------------
331 subroutine reada(rekord,lancuch,wartosc,default)
334 character*(*) :: rekord,lancuch
335 real(kind=8) :: wartosc,default
336 integer :: iread !,ilen
338 iread=index(rekord,lancuch)
343 iread=iread+ilen(lancuch)+1
344 read (rekord(iread:),*,err=10,end=10) wartosc
349 !-----------------------------------------------------------------------------
350 subroutine readi(rekord,lancuch,wartosc,default)
353 character*(*) :: rekord,lancuch
354 integer :: wartosc,default
355 integer :: iread !,ilen
357 iread=index(rekord,lancuch)
362 iread=iread+ilen(lancuch)+1
363 read (rekord(iread:),*,err=10,end=10) wartosc
368 !-----------------------------------------------------------------------------
369 subroutine multreadi(rekord,lancuch,tablica,dim,default)
373 integer :: tablica(dim),default
374 character*(*) :: rekord,lancuch
375 character(len=80) :: aux
376 integer :: iread !,ilen
381 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
382 if (iread.eq.0) return
383 iread=iread+ilen(lancuch)+1
384 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
386 end subroutine multreadi
387 !-----------------------------------------------------------------------------
388 subroutine multreada(rekord,lancuch,tablica,dim,default)
392 real(kind=8) :: tablica(dim),default
393 character*(*) :: rekord,lancuch
394 character(len=80) :: aux
395 integer :: iread !,ilen
400 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
401 if (iread.eq.0) return
402 iread=iread+ilen(lancuch)+1
403 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
405 end subroutine multreada
406 !-----------------------------------------------------------------------------
407 subroutine card_concat(card,to_upper)
409 ! dla UNRESA to_upper jest zawsze .true.
410 ! implicit real*8 (a-h,o-z)
411 ! include 'DIMENSIONS'
412 ! include 'COMMON.IOUNITS'
414 character(len=80) :: karta !,ucase
417 read (inp,'(a)') karta
418 if (to_upper) karta=ucase(karta)
420 do while (karta(80:80).eq.'&')
421 card=card(:ilen(card)+1)//karta(:79)
422 read (inp,'(a)') karta
423 if (to_upper) karta=ucase(karta)
425 card=card(:ilen(card)+1)//karta
427 end subroutine card_concat
428 !----------------------------------------------------------------------------
429 subroutine read_afminp
432 use control_data, only:out1file
434 character*320 afmcard
437 call card_concat(afmcard,.true.)
438 call readi(afmcard,"BEG",afmbeg,0)
439 call readi(afmcard,"END",afmend,0)
440 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
441 call reada(afmcard,"VEL",velAFMconst,0.0d0)
442 print *,'FORCE=' ,forceAFMconst
443 !------ NOW PROPERTIES FOR AFM
446 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
448 distafminit=dsqrt(distafminit)
449 print *,'initdist',distafminit
451 end subroutine read_afminp
452 !-----------------------------------------------------------------------------
453 subroutine read_dist_constr
456 use geometry, only: dist
460 ! implicit real*8 (a-h,o-z)
461 ! include 'DIMENSIONS'
465 ! include 'COMMON.SETUP'
466 ! include 'COMMON.CONTROL'
467 ! include 'COMMON.CHAIN'
468 ! include 'COMMON.IOUNITS'
469 ! include 'COMMON.SBRIDGE'
470 integer,dimension(2,100) :: ifrag_,ipair_
471 real(kind=8),dimension(100) :: wfrag_,wpair_
472 character(len=640) :: controlcard
475 integer :: i,k,j,ddjk,ii,jj,itemp
476 integer :: nfrag_,npair_,ndist_
477 real(kind=8) :: dist_cut
479 ! write (iout,*) "Calling read_dist_constr"
480 ! write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
482 if(.not.allocated(ihpb)) allocate(ihpb(maxdim))
483 if(.not.allocated(jhpb)) allocate(jhpb(maxdim))
484 if(.not.allocated(dhpb)) allocate(dhpb(maxdim))
485 if(.not.allocated(dhpb1)) allocate(dhpb1(maxdim))
486 if(.not.allocated(forcon)) allocate(forcon(maxdim))
487 if(.not.allocated(fordepth)) allocate(fordepth(maxdim))
488 if ((genconstr.gt.0).and.(constr_dist.eq.11)) then
489 call gen_dist_constr2
492 call card_concat(controlcard,.true.)
493 call readi(controlcard,"NFRAG",nfrag_,0)
494 call readi(controlcard,"NPAIR",npair_,0)
495 call readi(controlcard,"NDIST",ndist_,0)
496 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
497 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
498 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
499 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
500 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
501 ! write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
502 ! write (iout,*) "IFRAG"
504 ! write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
506 ! write (iout,*) "IPAIR"
508 ! write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
510 ! if(.not.allocated(ihpb)) allocate(ihpb(maxdim))
511 ! if(.not.allocated(jhpb)) allocate(jhpb(maxdim))
512 ! if(.not.allocated(dhpb)) allocate(dhpb(maxdim))
513 ! if(.not.allocated(forcon)) allocate(forcon(maxdim))
517 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
518 if (ifrag_(2,i).gt.nstart_sup+nsup-1) &
519 ifrag_(2,i)=nstart_sup+nsup-1
520 ! write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
522 if (wfrag_(i).gt.0.0d0) then
523 do j=ifrag_(1,i),ifrag_(2,i)-1
525 ! write (iout,*) "j",j," k",k
527 if (constr_dist.eq.1) then
532 forcon(nhpb)=wfrag_(i)
533 else if (constr_dist.eq.2) then
534 if (ddjk.le.dist_cut) then
539 forcon(nhpb)=wfrag_(i)
546 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
549 if (.not.out1file .or. me.eq.king) &
550 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",&
551 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
553 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",&
554 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
561 if (wpair_(i).gt.0.0d0) then
569 do j=ifrag_(1,ii),ifrag_(2,ii)
570 do k=ifrag_(1,jj),ifrag_(2,jj)
574 forcon(nhpb)=wpair_(i)
577 if (.not.out1file .or. me.eq.king) &
578 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
579 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
581 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
582 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
589 ! read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
590 ! if (forcon(nhpb+1).gt.0.0d0) then
592 ! dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
593 if (constr_dist.eq.11) then
594 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), &
595 ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
596 fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
599 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), &
600 ibecarb(i),forcon(nhpb+1)
602 if (forcon(nhpb+1).gt.0.0d0) then
604 if (ibecarb(i).gt.0) then
608 if (dhpb(nhpb).eq.0.0d0) &
609 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
613 if (.not.out1file .or. me.eq.king) &
614 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
615 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
617 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
618 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
624 end subroutine read_dist_constr
625 !-----------------------------------------------------------------------------
626 subroutine gen_dist_constr2
629 use geometry, only: dist
634 real(kind=8) :: distance
635 if (constr_dist.eq.11) then
636 do i=nstart_sup,nstart_sup+nsup-1
637 do j=i+2,nstart_sup+nsup-1
639 if (distance.le.15.0) then
641 ihpb(nhpb)=i+nstart_seq-nstart_sup
642 jhpb(nhpb)=j+nstart_seq-nstart_sup
643 forcon(nhpb)=sqrt(0.04*distance)
644 fordepth(nhpb)=sqrt(40.0/distance)
645 dhpb(nhpb)=distance-0.1d0
646 dhpb1(nhpb)=distance+0.1d0
649 if (.not.out1file .or. me.eq.king) &
650 write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ", &
651 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
653 write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ", &
654 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
660 do i=nstart_sup,nstart_sup+nsup-1
661 do j=i+2,nstart_sup+nsup-1
663 ihpb(nhpb)=i+nstart_seq-nstart_sup
664 jhpb(nhpb)=j+nstart_seq-nstart_sup
671 end subroutine gen_dist_constr2
673 !-----------------------------------------------------------------------------
685 !-----------------------------------------------------------------------------
686 subroutine copy_to_tmp(source)
688 ! include "DIMENSIONS"
689 ! include "COMMON.IOUNITS"
690 character*(*) :: source
691 character(len=256) :: tmpfile
695 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
696 inquire(file=tmpfile,exist=ex)
698 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),&
699 " to temporary directory..."
700 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
701 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
704 end subroutine copy_to_tmp
705 !-----------------------------------------------------------------------------
706 subroutine move_from_tmp(source)
708 ! include "DIMENSIONS"
709 ! include "COMMON.IOUNITS"
710 character*(*) :: source
713 write (*,*) "Moving ",source(:ilen(source)),&
714 " from temporary directory to working directory"
715 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
716 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
718 end subroutine move_from_tmp
719 !-----------------------------------------------------------------------------
721 !-----------------------------------------------------------------------------
722 ! $Date: 1994/10/12 17:24:21 $
725 logical function find_arg(ipos,line,errflag)
727 integer, parameter :: maxlen = 80
728 character(len=80) :: line
729 character(len=1) :: empty=' ',equal='='
732 ! This function returns .TRUE., if an argument follows keyword keywd; if so
733 ! IPOS will point to the first non-blank character of the argument. Returns
734 ! .FALSE., if no argument follows the keyword; in this case IPOS points
735 ! to the first non-blank character of the next keyword.
737 do while (line(ipos:ipos) .eq. empty .and. ipos.le.maxlen)
741 if (line(ipos:ipos).eq.equal) then
744 do while (line(ipos:ipos) .eq. empty .and. ipos.le.maxlen)
747 if (ipos.gt.maxlen) errflag=.true.
753 end function find_arg
754 !-----------------------------------------------------------------------------
755 logical function find_group(iunit,jout,key1)
757 character*(*) :: key1
758 character(len=80) :: karta !,ucase
759 integer :: iunit,jout
766 do while (index(ucase(karta),key1(1:ll)).eq.0.or.lcom(1,karta))
767 read (iunit,'(a)',end=10) karta
769 write (jout,'(2a)') '> ',karta(1:78)
772 10 find_group=.false.
774 end function find_group
775 !-----------------------------------------------------------------------------
776 logical function iblnk(charc)
777 character(len=1) :: charc
780 iblnk = (n.eq.9) .or. (n.eq.10) .or. (charc.eq.' ')
783 !-----------------------------------------------------------------------------
784 integer function ilen(string)
785 character*(*) :: string
789 1 if ( ilen .gt. 0 ) then
790 if ( iblnk( string(ilen:ilen) ) ) then
797 !-----------------------------------------------------------------------------
798 integer function in_keywd_set(nkey,ikey,narg,keywd,keywdset)
799 integer :: nkey,i,ikey,narg
800 character(len=16) :: keywd,keywdset(1:nkey,0:nkey)
801 ! character(len=16) :: ucase
804 if (ucase(keywd).eq.keywdset(i,ikey)) then
810 ! No match to the allowed set of keywords if this point is reached.
813 end function in_keywd_set
814 !-----------------------------------------------------------------------------
815 character function lcase(string)
816 integer :: i, k, idiff
817 character*(*) :: string
818 character(len=1) :: c
819 character(len=40) :: chtmp
825 if (string(k+1:) .ne. ' ') then
829 idiff = ichar('a') - ichar('A')
833 if (lge(c,'A') .and. lle(c,'Z')) then
834 lcase(i:i) = char(ichar(c) + idiff)
839 !-----------------------------------------------------------------------------
840 logical function lcom(ipos,karta)
841 character(len=80) :: karta
842 character :: koment(2) = (/'!','#'/)
847 if (karta(ipos:ipos).eq.koment(i)) lcom=.true.
851 !-----------------------------------------------------------------------------
852 logical function lower_case(ch)
854 lower_case=(ch.ge.'a' .and. ch.le.'z')
856 end function lower_case
857 !-----------------------------------------------------------------------------
858 subroutine mykey(line,keywd,ipos,blankline,errflag)
860 ! This subroutine seeks a non-empty substring keywd in the string LINE.
861 ! The substring begins with the first character different from blank and
862 ! "=" encountered right to the pointer IPOS (inclusively) and terminates
863 ! at the character left to the first blank or "=". When the subroutine is
864 ! exited, the pointer IPOS is moved to the position of the terminator in LINE.
865 ! The logical variable BLANKLINE is set at .TRUE., if LINE(IPOS:) contains
866 ! only separators or the maximum length of the data line (80) has been reached.
867 ! The logical variable ERRFLAG is set at .TRUE. if the string
868 ! consists only from a "=".
869 integer, parameter :: maxlen=80
870 character(len=1) :: empty=' ',equal='=',comma=','
871 character*(*) :: keywd
872 character(len=80) :: line
873 logical :: blankline,errflag !EL,lcom
874 integer :: ipos,istart,iend
876 do while (line(ipos:ipos).eq.empty .and. (ipos.le.maxlen))
879 if (ipos.gt.maxlen .or. lcom(ipos,line) ) then
880 ! At this point the rest of the input line turned out to contain only blanks
881 ! or to be commented out.
887 ! Checks whether the current char is a separator.
888 do while (line(ipos:ipos).ne.empty .and. line(ipos:ipos).ne.equal &
889 .and. line(ipos:ipos).ne.comma .and. ipos.le.maxlen)
893 ! Error flag set to .true., if the length of the keyword was found less than 1.
894 if (iend.lt.istart) then
898 keywd=line(istart:iend)
901 !-----------------------------------------------------------------------------
902 subroutine numstr(inum,numm)
903 character(len=10) :: huj='0123456789'
904 character*(*) :: numm
905 integer :: inumm,inum,inum1,inum2
910 numm(3:3)=huj(inum2+1:inum2+1)
914 numm(2:2)=huj(inum2+1:inum2+1)
918 numm(1:1)=huj(inum2+1:inum2+1)
920 end subroutine numstr
921 !-----------------------------------------------------------------------------
922 function ucase(string)
923 integer :: i, k, idiff
924 character(*) :: string
925 character(len=len(string)) :: ucase
926 character(len=1) :: c
927 character(len=40) :: chtmp
933 if (string(k+1:) .ne. ' ') then
937 idiff = ichar('a') - ichar('A')
941 if (lge(c,'a') .and. lle(c,'z')) then
942 ucase(i:i) = char(ichar(c) - idiff)
947 !-----------------------------------------------------------------------------
949 !-----------------------------------------------------------------------------
950 subroutine pdbout(etot,tytul,iunit)
952 use geometry_data, only: c,nres
957 ! implicit real*8 (a-h,o-z)
958 ! include 'DIMENSIONS'
959 ! include 'COMMON.CHAIN'
960 ! include 'COMMON.INTERACT'
961 ! include 'COMMON.NAMES'
962 ! include 'COMMON.IOUNITS'
963 ! include 'COMMON.HEADER'
964 ! include 'COMMON.SBRIDGE'
965 ! include 'COMMON.DISTFIT'
966 ! include 'COMMON.MD'
967 !el character(len=50) :: tytul
968 character*(*) :: tytul
969 character(len=1),dimension(10) :: chainid= (/'A','B','C','D','E','F','G','H','I','J'/)
970 integer,dimension(nres) :: ica !(maxres)
973 integer :: j,iti,itj,itk,itl,i,iatom,ichain,ires,iunit
978 if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
980 write (iunit,'(3a,1pe15.5)') 'REMARK ',tytul,' ENERGY ',etot
981 !model write (iunit,'(a5,i6)') 'MODEL',1
982 if (nhfrag.gt.0) then
984 iti=itype(hfrag(1,j),1)
985 itj=itype(hfrag(2,j),1)
987 write (iunit,'(a5,i5,1x,a1,i1,2x,a3,i7,2x,a3,i7,i3,t76,i5)') &
989 restyp(iti,1),hfrag(1,j)-1,&
990 restyp(itj,1),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j)
992 write (iunit,'(a5,i5,1x,a1,i2,1x,a3,i7,2x,a3,i7,i3)') &
994 restyp(iti,1),hfrag(1,j)-1,&
995 restyp(itj,1),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j)
1000 if (nbfrag.gt.0) then
1004 iti=itype(bfrag(1,j),1)
1005 itj=itype(bfrag(2,j)-1,1)
1007 write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3)') &
1009 restyp(iti,1),bfrag(1,j)-1,&
1010 restyp(itj,1),bfrag(2,j)-2,0
1012 if (bfrag(3,j).gt.bfrag(4,j)) then
1014 itk=itype(bfrag(3,j),1)
1015 itl=itype(bfrag(4,j)+1,1)
1017 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)') &
1019 restyp(itl,1),bfrag(4,j),&
1020 restyp(itk,1),bfrag(3,j)-1,-1,&
1021 "N",restyp(itk,1),bfrag(3,j)-1,&
1022 "O",restyp(iti,1),bfrag(1,j)-1
1026 itk=itype(bfrag(3,j),1)
1027 itl=itype(bfrag(4,j)-1,1)
1030 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)') &
1032 restyp(itk,1),bfrag(3,j)-1,&
1033 restyp(itl,1),bfrag(4,j)-2,1,&
1034 "N",restyp(itk,1),bfrag(3,j)-1,&
1035 "O",restyp(iti,1),bfrag(1,j)-1
1047 write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') &
1048 'SSBOND',i,'CYS',idssb(i)-nnt+1,&
1049 'CYS',jdssb(i)-nnt+1
1051 write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') &
1052 'SSBOND',i,'CYS',ihpb(i)-nnt+1-nres,&
1053 'CYS',jhpb(i)-nnt+1-nres
1064 if ((iti.eq.ntyp1).and.(iti1.eq.ntyp1)) cycle
1065 if (iti.eq.ntyp1) then
1068 write (iunit,'(a)') 'TER'
1073 write (iunit,10) iatom,restyp(iti,1),chainid(ichain),&
1074 ires,(c(j,i),j=1,3),vtot(i)
1077 write (iunit,20) iatom,restyp(iti,1),chainid(ichain),&
1078 ires,(c(j,nres+i),j=1,3),&
1083 write (iunit,'(a)') 'TER'
1085 if (itype(i,1).eq.ntyp1) cycle
1086 if (itype(i,1).eq.10 .and. itype(i+1,1).ne.ntyp1) then
1087 write (iunit,30) ica(i),ica(i+1)
1088 else if (itype(i,1).ne.10 .and. itype(i+1,1).ne.ntyp1) then
1089 write (iunit,30) ica(i),ica(i+1),ica(i)+1
1090 else if (itype(i,1).ne.10 .and. itype(i+1,1).eq.ntyp1) then
1091 write (iunit,30) ica(i),ica(i)+1
1094 if (itype(nct,1).ne.10) then
1095 write (iunit,30) ica(nct),ica(nct)+1
1099 write (iunit,30) ica(idssb(i))+1,ica(jdssb(i))+1
1101 write (iunit,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
1104 write (iunit,'(a6)') 'ENDMDL'
1105 10 FORMAT ('ATOM',I7,' CA ',A3,1X,A1,I4,4X,3F8.3,f15.3)
1106 20 FORMAT ('ATOM',I7,' CB ',A3,1X,A1,I4,4X,3F8.3,f15.3)
1107 30 FORMAT ('CONECT',8I5)
1109 end subroutine pdbout
1110 !-----------------------------------------------------------------------------
1111 subroutine MOL2out(etot,tytul)
1112 ! Prints the Cartesian coordinates of the alpha-carbons in the Tripos mol2
1114 use geometry_data, only: c
1117 ! implicit real*8 (a-h,o-z)
1118 ! include 'DIMENSIONS'
1119 ! include 'COMMON.CHAIN'
1120 ! include 'COMMON.INTERACT'
1121 ! include 'COMMON.NAMES'
1122 ! include 'COMMON.IOUNITS'
1123 ! include 'COMMON.HEADER'
1124 ! include 'COMMON.SBRIDGE'
1125 character(len=32) :: tytul,fd
1126 character(len=3) :: zahl
1127 character(len=6) :: res_num,pom !,ucase
1131 real(kind=8) :: etot
1135 #elif (defined CRAY)
1140 write (imol2,'(a)') '#'
1141 write (imol2,'(a)') &
1142 '# Creating user name: unres'
1143 write (imol2,'(2a)') '# Creation time: ',&
1145 write (imol2,'(/a)') '\@<TRIPOS>MOLECULE'
1146 write (imol2,'(a)') tytul
1147 write (imol2,'(5i5)') nct-nnt+1,nct-nnt+nss+1,nct-nnt+nss+1,0,0
1148 write (imol2,'(a)') 'SMALL'
1149 write (imol2,'(a)') 'USER_CHARGES'
1150 write (imol2,'(a)') '\@<TRIPOS>ATOM'
1152 write (zahl,'(i3)') i
1153 pom=ucase(restyp(itype(i,1),1))
1154 res_num = pom(:3)//zahl(2:)
1155 write (imol2,10) i,(c(j,i),j=1,3),i,res_num,0.0
1157 write (imol2,'(a)') '\@<TRIPOS>BOND'
1159 write (imol2,'(i5,2i6,i2)') i-nnt+1,i-nnt+1,i-nnt+2,1
1162 write (imol2,'(i5,2i6,i2)') nct-nnt+i,ihpb(i),jhpb(i),1
1164 write (imol2,'(a)') '\@<TRIPOS>SUBSTRUCTURE'
1166 write (zahl,'(i3)') i
1167 pom = ucase(restyp(itype(i,1),1))
1168 res_num = pom(:3)//zahl(2:)
1169 write (imol2,30) i-nnt+1,res_num,i-nnt+1,0
1171 10 FORMAT (I7,' CA ',3F10.4,' C.3',I8,1X,A,F11.4,' ****')
1172 30 FORMAT (I7,1x,A,I14,' RESIDUE',I13,' **** ****')
1174 end subroutine MOL2out
1175 !-----------------------------------------------------------------------------
1179 use energy_data, only: itype
1181 ! implicit real*8 (a-h,o-z)
1182 ! include 'DIMENSIONS'
1183 ! include 'COMMON.IOUNITS'
1184 ! include 'COMMON.CHAIN'
1185 ! include 'COMMON.VAR'
1186 ! include 'COMMON.LOCAL'
1187 ! include 'COMMON.INTERACT'
1188 ! include 'COMMON.NAMES'
1189 ! include 'COMMON.GEO'
1190 ! include 'COMMON.TORSION'
1194 write (iout,'(/a)') 'Geometry of the virtual chain.'
1195 write (iout,'(7a)') ' Res ',' d',' Theta',&
1196 ' Phi',' Dsc',' Alpha',' Omega'
1199 ! print *,vbld(i),"vbld(i)",i
1200 write (iout,'(a3,i4,6f10.3)') restyp(iti,1),i,vbld(i),&
1201 rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),rad2deg*alph(i),&
1205 end subroutine intout
1206 !-----------------------------------------------------------------------------
1208 subroutine briefout(it,ener,free)!,plik)
1210 subroutine briefout(it,ener)
1215 ! implicit real*8 (a-h,o-z)
1216 ! include 'DIMENSIONS'
1217 ! include 'COMMON.IOUNITS'
1218 ! include 'COMMON.CHAIN'
1219 ! include 'COMMON.VAR'
1220 ! include 'COMMON.LOCAL'
1221 ! include 'COMMON.INTERACT'
1222 ! include 'COMMON.NAMES'
1223 ! include 'COMMON.GEO'
1224 ! include 'COMMON.SBRIDGE'
1225 ! print '(a,i5)',intname,igeom
1228 real(kind=8) :: ener,free
1229 ! character(len=80) :: plik
1232 #if defined(AIX) || defined(PGI)
1233 open (igeom,file=intname,position='append')
1235 open (igeom,file=intname,access='append')
1244 WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,NSS)
1246 WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,8)
1248 WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,NSS)
1250 WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,9)
1252 WRITE (igeom,190) (IHPB(I),JHPB(I),I=10,NSS)
1254 ! IF (nvar.gt.nphi) WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
1255 WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
1256 WRITE (igeom,200) (RAD2DEG*PHI(I),I=4,NRES)
1257 ! if (nvar.gt.nphi+ntheta) then
1258 write (igeom,200) (rad2deg*alph(i),i=2,nres-1)
1259 write (igeom,200) (rad2deg*omeg(i),i=2,nres-1)
1262 180 format (I5,F12.3,I2,9(1X,2I3))
1263 190 format (3X,11(1X,2I3))
1266 end subroutine briefout
1267 !-----------------------------------------------------------------------------
1269 subroutine fdate(fd)
1270 character(len=32) :: fd
1273 end subroutine fdate
1275 !-----------------------------------------------------------------------------
1277 real(kind=8) function gyrate(jcon)
1279 real(kind=8) function gyrate()
1282 use geometry_data, only: c
1285 ! implicit real*8 (a-h,o-z)
1286 ! include 'DIMENSIONS'
1287 ! include 'COMMON.INTERACT'
1288 ! include 'COMMON.CHAIN'
1290 real(kind=8),dimension(3) :: cen
1300 cen(j)=cen(j)+c(j,i)
1304 cen(j)=cen(j)/dble(nct-nnt+1)
1309 rg = rg + (c(j,i)-cen(j))**2
1313 gyrate = dsqrt(rg/dble(nct-nnt+1))
1315 gyrate = sqrt(rg/dble(nct-nnt+1))
1320 !-----------------------------------------------------------------------------
1322 subroutine reads(rekord,lancuch,wartosc,default)
1324 character*(*) :: rekord,lancuch,wartosc,default
1325 character(len=80) :: aux
1326 integer :: lenlan,lenrec,iread,ireade
1330 lenlan=ilen(lancuch)
1332 iread=index(rekord,lancuch(:lenlan)//"=")
1333 ! print *,"rekord",rekord," lancuch",lancuch
1334 ! print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
1335 if (iread.eq.0) then
1339 iread=iread+lenlan+1
1340 ! print *,"iread",iread
1341 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
1342 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
1344 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
1346 ! print *,"iread",iread
1347 if (iread.gt.lenrec) then
1352 ! print *,"ireade",ireade
1353 do while (ireade.lt.lenrec .and. &
1354 .not.iblnk(rekord(ireade:ireade)))
1357 wartosc=rekord(iread:ireade)
1359 end subroutine reads
1361 !-----------------------------------------------------------------------------
1363 !-----------------------------------------------------------------------------
1364 subroutine permut(isym)
1366 use geometry_data, only: tabperm
1367 ! implicit real*8 (a-h,o-z)
1368 ! include 'DIMENSIONS'
1369 ! include 'COMMON.LOCAL'
1370 ! include 'COMMON.VAR'
1371 ! include 'COMMON.CHAIN'
1372 ! include 'COMMON.INTERACT'
1373 ! include 'COMMON.IOUNITS'
1374 ! include 'COMMON.GEO'
1375 ! include 'COMMON.NAMES'
1376 ! include 'COMMON.CONTROL'
1381 integer,dimension(isym) :: a
1382 ! parameter(n=symetr)
1395 10 print *,(a(i),i=1,n)
1399 ! write (iout,*) "tututu", kkk
1401 if(nextp(n,a)) go to 10
1403 end subroutine permut
1404 !-----------------------------------------------------------------------------
1405 logical function nextp(n,a)
1407 integer :: n,i,j,k,t
1409 integer,dimension(n) :: a
1411 10 if(a(i).lt.a(i+1)) go to 20
1428 if(a(j).lt.a(i)) go to 40
1435 !-----------------------------------------------------------------------------
1436 !-----------------------------------------------------------------------------