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
1062 iti=itype(i,molnum(i))
1064 if (molnum(i+1).eq.0) then
1065 iti1=ntyp1_molec(molnum(i))
1067 iti1=itype(i+1,molnum(i+1))
1069 if ((iti.eq.ntyp1_molec(molnum(i))).and.(iti1.eq.ntyp1_molec(molnum(i)))) cycle
1071 if (iti.eq.ntyp1_molec(molnum(i)).and.(molnum(i+1).eq.5)) cycle
1073 if (iti.eq.ntyp1_molec(molnum(i))) then
1076 write (iunit,'(a)') 'TER'
1081 if (molnum(i).eq.1) then
1083 write (iunit,10) iatom,restyp(iti,molnum(i)),chainid(ichain),&
1084 ires,(c(j,i),j=1,3),vtot(i)
1085 elseif(molnum(i).eq.2) then
1086 if (istype(i).eq.0) istype(i)=1
1087 write (iunit,40) iatom,sugartyp(istype(i)),restyp(iti,2), &
1088 chainid(ichain),ires,(c(j,i),j=1,3),vtot(i)
1090 write (iunit,60) iatom,restyp(iti,molnum(i)),chainid(ichain),&
1091 ires,(c(j,i),j=1,3),vtot(i)
1093 if ((iti.ne.10).and.(molnum(i).ne.5)) then
1095 if (molnum(i).eq.1) then
1096 write (iunit,20) iatom,restyp(iti,1),chainid(ichain),&
1097 ires,(c(j,nres+i),j=1,3),&
1099 else if (molnum(i).eq.2) then
1100 if (istype(i).eq.0) istype(i)=1
1101 write (iunit,50) iatom,sugartyp(istype(i)),restyp(iti,2), &
1102 chainid(ichain),ires,(c(j,nres+i),j=1,3),vtot(i+nres)
1108 write (iunit,'(a)') 'TER'
1110 if (itype(i,1).eq.ntyp1) cycle
1111 if ((itype(i,1).eq.10 .and. itype(i+1,1).ne.ntyp1).or.(molnum(i).eq.5)) then
1112 write (iunit,30) ica(i),ica(i+1)
1113 else if (itype(i,1).ne.10 .and. itype(i+1,1).ne.ntyp1) then
1114 write (iunit,30) ica(i),ica(i+1),ica(i)+1
1115 else if (itype(i,1).ne.10 .and. itype(i+1,1).eq.ntyp1) then
1116 write (iunit,30) ica(i),ica(i)+1
1119 if ((itype(nct,1).ne.10).and.(molnum(i).ne.5)) then
1120 write (iunit,30) ica(nct),ica(nct)+1
1124 write (iunit,30) ica(idssb(i))+1,ica(jdssb(i))+1
1126 write (iunit,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
1129 write (iunit,'(a6)') 'ENDMDL'
1130 10 FORMAT ('ATOM',I7,' CA ',A3,1X,A1,I4,4X,3F8.3,f15.3)
1132 20 FORMAT ('ATOM',I7,' CB ',A3,1X,A1,I4,4X,3F8.3,f15.3)
1133 40 FORMAT ("ATOM",I7," C5' ",1X,2A1,1X,A1,I4,4X,3F8.3,f15.3)
1134 50 FORMAT ("ATOM",I7," C1' ",1X,2A1,1X,A1,I4,4X,3F8.3,f15.3)
1136 30 FORMAT ('CONECT',8I5)
1137 60 FORMAT ('HETATM',I5,' CA ',A3,1X,A1,I4,4X,3F8.3,f15.3)
1140 end subroutine pdbout
1141 !-----------------------------------------------------------------------------
1142 subroutine MOL2out(etot,tytul)
1143 ! Prints the Cartesian coordinates of the alpha-carbons in the Tripos mol2
1145 use geometry_data, only: c
1148 ! implicit real*8 (a-h,o-z)
1149 ! include 'DIMENSIONS'
1150 ! include 'COMMON.CHAIN'
1151 ! include 'COMMON.INTERACT'
1152 ! include 'COMMON.NAMES'
1153 ! include 'COMMON.IOUNITS'
1154 ! include 'COMMON.HEADER'
1155 ! include 'COMMON.SBRIDGE'
1156 character(len=32) :: tytul,fd
1157 character(len=3) :: zahl
1158 character(len=6) :: res_num,pom !,ucase
1162 real(kind=8) :: etot
1166 #elif (defined CRAY)
1171 write (imol2,'(a)') '#'
1172 write (imol2,'(a)') &
1173 '# Creating user name: unres'
1174 write (imol2,'(2a)') '# Creation time: ',&
1176 write (imol2,'(/a)') '\@<TRIPOS>MOLECULE'
1177 write (imol2,'(a)') tytul
1178 write (imol2,'(5i5)') nct-nnt+1,nct-nnt+nss+1,nct-nnt+nss+1,0,0
1179 write (imol2,'(a)') 'SMALL'
1180 write (imol2,'(a)') 'USER_CHARGES'
1181 write (imol2,'(a)') '\@<TRIPOS>ATOM'
1183 write (zahl,'(i3)') i
1184 pom=ucase(restyp(itype(i,1),1))
1185 res_num = pom(:3)//zahl(2:)
1186 write (imol2,10) i,(c(j,i),j=1,3),i,res_num,0.0
1188 write (imol2,'(a)') '\@<TRIPOS>BOND'
1190 write (imol2,'(i5,2i6,i2)') i-nnt+1,i-nnt+1,i-nnt+2,1
1193 write (imol2,'(i5,2i6,i2)') nct-nnt+i,ihpb(i),jhpb(i),1
1195 write (imol2,'(a)') '\@<TRIPOS>SUBSTRUCTURE'
1197 write (zahl,'(i3)') i
1198 pom = ucase(restyp(itype(i,1),1))
1199 res_num = pom(:3)//zahl(2:)
1200 write (imol2,30) i-nnt+1,res_num,i-nnt+1,0
1202 10 FORMAT (I7,' CA ',3F10.4,' C.3',I8,1X,A,F11.4,' ****')
1203 30 FORMAT (I7,1x,A,I14,' RESIDUE',I13,' **** ****')
1205 end subroutine MOL2out
1206 !-----------------------------------------------------------------------------
1210 use energy_data, only: itype
1212 ! implicit real*8 (a-h,o-z)
1213 ! include 'DIMENSIONS'
1214 ! include 'COMMON.IOUNITS'
1215 ! include 'COMMON.CHAIN'
1216 ! include 'COMMON.VAR'
1217 ! include 'COMMON.LOCAL'
1218 ! include 'COMMON.INTERACT'
1219 ! include 'COMMON.NAMES'
1220 ! include 'COMMON.GEO'
1221 ! include 'COMMON.TORSION'
1225 write (iout,'(/a)') 'Geometry of the virtual chain.'
1226 write (iout,'(7a)') ' Res ',' d',' Theta',&
1227 ' Phi',' Dsc',' Alpha',' Omega'
1230 ! print *,vbld(i),"vbld(i)",i
1231 write (iout,'(a3,i4,6f10.3)') restyp(iti,1),i,vbld(i),&
1232 rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),rad2deg*alph(i),&
1236 end subroutine intout
1237 !-----------------------------------------------------------------------------
1239 subroutine briefout(it,ener,free)!,plik)
1241 subroutine briefout(it,ener)
1246 ! implicit real*8 (a-h,o-z)
1247 ! include 'DIMENSIONS'
1248 ! include 'COMMON.IOUNITS'
1249 ! include 'COMMON.CHAIN'
1250 ! include 'COMMON.VAR'
1251 ! include 'COMMON.LOCAL'
1252 ! include 'COMMON.INTERACT'
1253 ! include 'COMMON.NAMES'
1254 ! include 'COMMON.GEO'
1255 ! include 'COMMON.SBRIDGE'
1256 ! print '(a,i5)',intname,igeom
1259 real(kind=8) :: ener,free
1260 ! character(len=80) :: plik
1263 #if defined(AIX) || defined(PGI)
1264 open (igeom,file=intname,position='append')
1266 open (igeom,file=intname,access='append')
1275 WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,NSS)
1277 WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,8)
1279 WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,NSS)
1281 WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,9)
1283 WRITE (igeom,190) (IHPB(I),JHPB(I),I=10,NSS)
1285 ! IF (nvar.gt.nphi) WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
1286 WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
1287 WRITE (igeom,200) (RAD2DEG*PHI(I),I=4,NRES)
1288 ! if (nvar.gt.nphi+ntheta) then
1289 write (igeom,200) (rad2deg*alph(i),i=2,nres-1)
1290 write (igeom,200) (rad2deg*omeg(i),i=2,nres-1)
1293 180 format (I5,F12.3,I2,9(1X,2I3))
1294 190 format (3X,11(1X,2I3))
1297 end subroutine briefout
1298 !-----------------------------------------------------------------------------
1300 subroutine fdate(fd)
1301 character(len=32) :: fd
1304 end subroutine fdate
1306 !-----------------------------------------------------------------------------
1308 real(kind=8) function gyrate(jcon)
1310 real(kind=8) function gyrate()
1313 use geometry_data, only: c
1316 ! implicit real*8 (a-h,o-z)
1317 ! include 'DIMENSIONS'
1318 ! include 'COMMON.INTERACT'
1319 ! include 'COMMON.CHAIN'
1321 real(kind=8),dimension(3) :: cen
1331 cen(j)=cen(j)+c(j,i)
1335 cen(j)=cen(j)/dble(nct-nnt+1)
1340 rg = rg + (c(j,i)-cen(j))**2
1344 gyrate = dsqrt(rg/dble(nct-nnt+1))
1346 gyrate = sqrt(rg/dble(nct-nnt+1))
1351 !-----------------------------------------------------------------------------
1353 subroutine reads(rekord,lancuch,wartosc,default)
1355 character*(*) :: rekord,lancuch,wartosc,default
1356 character(len=80) :: aux
1357 integer :: lenlan,lenrec,iread,ireade
1361 lenlan=ilen(lancuch)
1363 iread=index(rekord,lancuch(:lenlan)//"=")
1364 ! print *,"rekord",rekord," lancuch",lancuch
1365 ! print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
1366 if (iread.eq.0) then
1370 iread=iread+lenlan+1
1371 ! print *,"iread",iread
1372 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
1373 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
1375 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
1377 ! print *,"iread",iread
1378 if (iread.gt.lenrec) then
1383 ! print *,"ireade",ireade
1384 do while (ireade.lt.lenrec .and. &
1385 .not.iblnk(rekord(ireade:ireade)))
1388 wartosc=rekord(iread:ireade)
1390 end subroutine reads
1392 !-----------------------------------------------------------------------------
1394 !-----------------------------------------------------------------------------
1395 subroutine permut(isym)
1397 use geometry_data, only: tabperm
1398 ! implicit real*8 (a-h,o-z)
1399 ! include 'DIMENSIONS'
1400 ! include 'COMMON.LOCAL'
1401 ! include 'COMMON.VAR'
1402 ! include 'COMMON.CHAIN'
1403 ! include 'COMMON.INTERACT'
1404 ! include 'COMMON.IOUNITS'
1405 ! include 'COMMON.GEO'
1406 ! include 'COMMON.NAMES'
1407 ! include 'COMMON.CONTROL'
1412 integer,dimension(isym) :: a
1413 ! parameter(n=symetr)
1426 10 print *,(a(i),i=1,n)
1430 ! write (iout,*) "tututu", kkk
1432 if(nextp(n,a)) go to 10
1434 end subroutine permut
1435 !-----------------------------------------------------------------------------
1436 logical function nextp(n,a)
1438 integer :: n,i,j,k,t
1440 integer,dimension(n) :: a
1442 10 if(a(i).lt.a(i+1)) go to 20
1459 if(a(j).lt.a(i)) go to 40
1466 !-----------------------------------------------------------------------------
1467 !-----------------------------------------------------------------------------