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) then
1076 write (iunit,'(a)') 'TER'
1081 if (molnum(i).ne.5) then
1082 write (iunit,10) iatom,restyp(iti,molnum(i)),chainid(ichain),&
1083 ires,(c(j,i),j=1,3),vtot(i)
1085 write (iunit,60) iatom,restyp(iti,molnum(i)),chainid(ichain),&
1086 ires,(c(j,i),j=1,3),vtot(i)
1088 if ((iti.ne.10).and.(molnum(i).ne.5)) then
1090 write (iunit,20) iatom,restyp(iti,1),chainid(ichain),&
1091 ires,(c(j,nres+i),j=1,3),&
1096 write (iunit,'(a)') 'TER'
1098 if (itype(i,1).eq.ntyp1) cycle
1099 if ((itype(i,1).eq.10 .and. itype(i+1,1).ne.ntyp1).or.(molnum(i).eq.5)) then
1100 write (iunit,30) ica(i),ica(i+1)
1101 else if (itype(i,1).ne.10 .and. itype(i+1,1).ne.ntyp1) then
1102 write (iunit,30) ica(i),ica(i+1),ica(i)+1
1103 else if (itype(i,1).ne.10 .and. itype(i+1,1).eq.ntyp1) then
1104 write (iunit,30) ica(i),ica(i)+1
1107 if ((itype(nct,1).ne.10).and.(molnum(i).ne.5)) then
1108 write (iunit,30) ica(nct),ica(nct)+1
1112 write (iunit,30) ica(idssb(i))+1,ica(jdssb(i))+1
1114 write (iunit,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
1117 write (iunit,'(a6)') 'ENDMDL'
1118 10 FORMAT ('ATOM',I7,' CA ',A3,1X,A1,I4,4X,3F8.3,f15.3)
1120 20 FORMAT ('ATOM',I7,' CB ',A3,1X,A1,I4,4X,3F8.3,f15.3)
1121 30 FORMAT ('CONECT',8I5)
1122 60 FORMAT ('HETATM',I5,' CA ',A3,1X,A1,I4,4X,3F8.3,f15.3)
1125 end subroutine pdbout
1126 !-----------------------------------------------------------------------------
1127 subroutine MOL2out(etot,tytul)
1128 ! Prints the Cartesian coordinates of the alpha-carbons in the Tripos mol2
1130 use geometry_data, only: c
1133 ! implicit real*8 (a-h,o-z)
1134 ! include 'DIMENSIONS'
1135 ! include 'COMMON.CHAIN'
1136 ! include 'COMMON.INTERACT'
1137 ! include 'COMMON.NAMES'
1138 ! include 'COMMON.IOUNITS'
1139 ! include 'COMMON.HEADER'
1140 ! include 'COMMON.SBRIDGE'
1141 character(len=32) :: tytul,fd
1142 character(len=3) :: zahl
1143 character(len=6) :: res_num,pom !,ucase
1147 real(kind=8) :: etot
1151 #elif (defined CRAY)
1156 write (imol2,'(a)') '#'
1157 write (imol2,'(a)') &
1158 '# Creating user name: unres'
1159 write (imol2,'(2a)') '# Creation time: ',&
1161 write (imol2,'(/a)') '\@<TRIPOS>MOLECULE'
1162 write (imol2,'(a)') tytul
1163 write (imol2,'(5i5)') nct-nnt+1,nct-nnt+nss+1,nct-nnt+nss+1,0,0
1164 write (imol2,'(a)') 'SMALL'
1165 write (imol2,'(a)') 'USER_CHARGES'
1166 write (imol2,'(a)') '\@<TRIPOS>ATOM'
1168 write (zahl,'(i3)') i
1169 pom=ucase(restyp(itype(i,1),1))
1170 res_num = pom(:3)//zahl(2:)
1171 write (imol2,10) i,(c(j,i),j=1,3),i,res_num,0.0
1173 write (imol2,'(a)') '\@<TRIPOS>BOND'
1175 write (imol2,'(i5,2i6,i2)') i-nnt+1,i-nnt+1,i-nnt+2,1
1178 write (imol2,'(i5,2i6,i2)') nct-nnt+i,ihpb(i),jhpb(i),1
1180 write (imol2,'(a)') '\@<TRIPOS>SUBSTRUCTURE'
1182 write (zahl,'(i3)') i
1183 pom = ucase(restyp(itype(i,1),1))
1184 res_num = pom(:3)//zahl(2:)
1185 write (imol2,30) i-nnt+1,res_num,i-nnt+1,0
1187 10 FORMAT (I7,' CA ',3F10.4,' C.3',I8,1X,A,F11.4,' ****')
1188 30 FORMAT (I7,1x,A,I14,' RESIDUE',I13,' **** ****')
1190 end subroutine MOL2out
1191 !-----------------------------------------------------------------------------
1195 use energy_data, only: itype
1197 ! implicit real*8 (a-h,o-z)
1198 ! include 'DIMENSIONS'
1199 ! include 'COMMON.IOUNITS'
1200 ! include 'COMMON.CHAIN'
1201 ! include 'COMMON.VAR'
1202 ! include 'COMMON.LOCAL'
1203 ! include 'COMMON.INTERACT'
1204 ! include 'COMMON.NAMES'
1205 ! include 'COMMON.GEO'
1206 ! include 'COMMON.TORSION'
1210 write (iout,'(/a)') 'Geometry of the virtual chain.'
1211 write (iout,'(7a)') ' Res ',' d',' Theta',&
1212 ' Phi',' Dsc',' Alpha',' Omega'
1215 ! print *,vbld(i),"vbld(i)",i
1216 write (iout,'(a3,i4,6f10.3)') restyp(iti,1),i,vbld(i),&
1217 rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),rad2deg*alph(i),&
1221 end subroutine intout
1222 !-----------------------------------------------------------------------------
1224 subroutine briefout(it,ener,free)!,plik)
1226 subroutine briefout(it,ener)
1231 ! implicit real*8 (a-h,o-z)
1232 ! include 'DIMENSIONS'
1233 ! include 'COMMON.IOUNITS'
1234 ! include 'COMMON.CHAIN'
1235 ! include 'COMMON.VAR'
1236 ! include 'COMMON.LOCAL'
1237 ! include 'COMMON.INTERACT'
1238 ! include 'COMMON.NAMES'
1239 ! include 'COMMON.GEO'
1240 ! include 'COMMON.SBRIDGE'
1241 ! print '(a,i5)',intname,igeom
1244 real(kind=8) :: ener,free
1245 ! character(len=80) :: plik
1248 #if defined(AIX) || defined(PGI)
1249 open (igeom,file=intname,position='append')
1251 open (igeom,file=intname,access='append')
1260 WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,NSS)
1262 WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,8)
1264 WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,NSS)
1266 WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,9)
1268 WRITE (igeom,190) (IHPB(I),JHPB(I),I=10,NSS)
1270 ! IF (nvar.gt.nphi) WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
1271 WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
1272 WRITE (igeom,200) (RAD2DEG*PHI(I),I=4,NRES)
1273 ! if (nvar.gt.nphi+ntheta) then
1274 write (igeom,200) (rad2deg*alph(i),i=2,nres-1)
1275 write (igeom,200) (rad2deg*omeg(i),i=2,nres-1)
1278 180 format (I5,F12.3,I2,9(1X,2I3))
1279 190 format (3X,11(1X,2I3))
1282 end subroutine briefout
1283 !-----------------------------------------------------------------------------
1285 subroutine fdate(fd)
1286 character(len=32) :: fd
1289 end subroutine fdate
1291 !-----------------------------------------------------------------------------
1293 real(kind=8) function gyrate(jcon)
1295 real(kind=8) function gyrate()
1298 use geometry_data, only: c
1301 ! implicit real*8 (a-h,o-z)
1302 ! include 'DIMENSIONS'
1303 ! include 'COMMON.INTERACT'
1304 ! include 'COMMON.CHAIN'
1306 real(kind=8),dimension(3) :: cen
1316 cen(j)=cen(j)+c(j,i)
1320 cen(j)=cen(j)/dble(nct-nnt+1)
1325 rg = rg + (c(j,i)-cen(j))**2
1329 gyrate = dsqrt(rg/dble(nct-nnt+1))
1331 gyrate = sqrt(rg/dble(nct-nnt+1))
1336 !-----------------------------------------------------------------------------
1338 subroutine reads(rekord,lancuch,wartosc,default)
1340 character*(*) :: rekord,lancuch,wartosc,default
1341 character(len=80) :: aux
1342 integer :: lenlan,lenrec,iread,ireade
1346 lenlan=ilen(lancuch)
1348 iread=index(rekord,lancuch(:lenlan)//"=")
1349 ! print *,"rekord",rekord," lancuch",lancuch
1350 ! print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
1351 if (iread.eq.0) then
1355 iread=iread+lenlan+1
1356 ! print *,"iread",iread
1357 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
1358 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
1360 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
1362 ! print *,"iread",iread
1363 if (iread.gt.lenrec) then
1368 ! print *,"ireade",ireade
1369 do while (ireade.lt.lenrec .and. &
1370 .not.iblnk(rekord(ireade:ireade)))
1373 wartosc=rekord(iread:ireade)
1375 end subroutine reads
1377 !-----------------------------------------------------------------------------
1379 !-----------------------------------------------------------------------------
1380 subroutine permut(isym)
1382 use geometry_data, only: tabperm
1383 ! implicit real*8 (a-h,o-z)
1384 ! include 'DIMENSIONS'
1385 ! include 'COMMON.LOCAL'
1386 ! include 'COMMON.VAR'
1387 ! include 'COMMON.CHAIN'
1388 ! include 'COMMON.INTERACT'
1389 ! include 'COMMON.IOUNITS'
1390 ! include 'COMMON.GEO'
1391 ! include 'COMMON.NAMES'
1392 ! include 'COMMON.CONTROL'
1397 integer,dimension(isym) :: a
1398 ! parameter(n=symetr)
1411 10 print *,(a(i),i=1,n)
1415 ! write (iout,*) "tututu", kkk
1417 if(nextp(n,a)) go to 10
1419 end subroutine permut
1420 !-----------------------------------------------------------------------------
1421 logical function nextp(n,a)
1423 integer :: n,i,j,k,t
1425 integer,dimension(n) :: a
1427 10 if(a(i).lt.a(i+1)) go to 20
1444 if(a(j).lt.a(i)) go to 40
1451 !-----------------------------------------------------------------------------
1452 !-----------------------------------------------------------------------------