2 !-----------------------------------------------------------------------
6 !-----------------------------------------------------------------------------
7 ! Max. number of AA residues
8 integer,parameter :: maxres=6000!1200
9 ! Appr. max. number of interaction sites
10 integer,parameter :: maxres2=2*maxres
11 ! parameter (maxres6=6*maxres)
12 ! parameter (mmaxres2=(maxres2*(maxres2+1)/2))
13 !-----------------------------------------------------------------------------
14 ! Max. number of S-S bridges
15 ! integer,parameter :: maxss=20
16 !-----------------------------------------------------------------------------
17 ! Max. number of derivatives of virtual-bond and side-chain vectors in theta
19 ! integer,parameter :: maxdim=(maxres-1)*(maxres-2)/2
20 !-----------------------------------------------------------------------------
23 !-----------------------------------------------------------------------------
25 !-----------------------------------------------------------------------------
27 !-----------------------------------------------------------------------------
28 subroutine read_bridge
29 ! Read information about disulfide bridges.
30 use geometry_data, only: nres
32 use control_data, only:out1file
34 ! implicit real*8 (a-h,o-z)
35 ! include 'DIMENSIONS'
39 ! include 'COMMON.IOUNITS'
40 ! include 'COMMON.GEO'
41 ! include 'COMMON.VAR'
42 ! include 'COMMON.INTERACT'
43 ! include 'COMMON.LOCAL'
44 ! include 'COMMON.NAMES'
45 ! include 'COMMON.CHAIN'
46 ! include 'COMMON.FFIELD'
47 ! include 'COMMON.SBRIDGE'
48 ! include 'COMMON.HEADER'
49 ! include 'COMMON.CONTROL'
50 ! include 'COMMON.DBASE'
51 ! include 'COMMON.THREAD'
52 ! include 'COMMON.TIME1'
53 ! include 'COMMON.SETUP'
57 ! Read bridging residues.
61 read (inp,*) (iss(i),i=1,ns)
65 if(me.eq.king.or..not.out1file) &
66 write (iout,*) 'ns=',ns
68 write(iout,*) ' iss:',(iss(i),i=1,ns)
69 ! Check whether the specified bridging residues are cystines.
71 write(iout,*) i,iss(i)
72 if (itype(iss(i),1).ne.1) then
73 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)') &
74 'Do you REALLY think that the residue ',&
75 restyp(itype(iss(i),1),1),i,&
76 ' can form a disulfide bridge?!!!'
77 write (*,'(2a,i3,a)') &
78 'Do you REALLY think that the residue ',&
79 restyp(itype(iss(i),1),1),i,&
80 ' can form a disulfide bridge?!!!'
82 call MPI_Finalize(MPI_COMM_WORLD,ierror)
88 if(.not.allocated(ihpb)) allocate(ihpb(ns))
89 if(.not.allocated(jhpb)) allocate(jhpb(ns))
91 ! Read preformed bridges.
95 if(.not.allocated(ihpb)) allocate(ihpb(nss))
96 if(.not.allocated(jhpb)) allocate(jhpb(nss))
97 read (inp,*) (ihpb(i),jhpb(i),i=1,nss)
99 if(.not.allocated(dhpb)) allocate(dhpb(nss))
100 if(.not.allocated(forcon)) allocate(forcon(nss))!(maxdim) !el maxdim=(maxres-1)*(maxres-2)/2
101 if(.not.allocated(dhpb1)) allocate(dhpb1(nss))
102 if(.not.allocated(ibecarb)) allocate(ibecarb(nss))
103 ! el Initialize the bridge array
108 !--------------------
110 write(iout,*)'nss=',nss !el,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
112 write(iout,*)'ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
114 ! Check if the residues involved in bridges are in the specified list of
118 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j) &
119 .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
120 write (iout,'(a,i3,a)') 'Disulfide pair',i,&
121 ' contains residues present in other pairs.'
122 write (*,'(a,i3,a)') 'Disulfide pair',i,&
123 ' contains residues present in other pairs.'
125 call MPI_Finalize(MPI_COMM_WORLD,ierror)
131 if (ihpb(i).eq.iss(j)) goto 10
133 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
136 if (jhpb(i).eq.iss(j)) goto 20
138 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
149 ! write(iout,*) "end read_bridge"
151 end subroutine read_bridge
152 !-----------------------------------------------------------------------------
153 subroutine read_x(kanal,*)
157 use geometry, only:int_from_cart1
158 ! implicit real*8 (a-h,o-z)
159 ! include 'DIMENSIONS'
160 ! include 'COMMON.GEO'
161 ! include 'COMMON.VAR'
162 ! include 'COMMON.CHAIN'
163 ! include 'COMMON.IOUNITS'
164 ! include 'COMMON.CONTROL'
165 ! include 'COMMON.LOCAL'
166 ! include 'COMMON.INTERACT'
167 ! Read coordinates from input
170 integer :: l,k,j,i,kanal
172 read(kanal,'(8f10.5)',end=10,err=10) &
173 ((c(l,k),l=1,3),k=1,nres),&
174 ((c(l,k+nres),l=1,3),k=nnt,nct)
177 c(j,2*nres)=c(j,nres)
179 call int_from_cart1(.false.)
182 dc(j,i)=c(j,i+1)-c(j,i)
183 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
187 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
189 dc(j,i+nres)=c(j,i+nres)-c(j,i)
190 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
197 end subroutine read_x
198 !-----------------------------------------------------------------------------
199 subroutine read_threadbase
204 ! implicit real*8 (a-h,o-z)
205 ! include 'DIMENSIONS'
206 ! include 'COMMON.IOUNITS'
207 ! include 'COMMON.GEO'
208 ! include 'COMMON.VAR'
209 ! include 'COMMON.INTERACT'
210 ! include 'COMMON.LOCAL'
211 ! include 'COMMON.NAMES'
212 ! include 'COMMON.CHAIN'
213 ! include 'COMMON.FFIELD'
214 ! include 'COMMON.SBRIDGE'
215 ! include 'COMMON.HEADER'
216 ! include 'COMMON.CONTROL'
217 ! include 'COMMON.DBASE'
218 ! include 'COMMON.THREAD'
219 ! include 'COMMON.TIME1'
224 ! Read pattern database for threading.
226 allocate(cart_base(3,maxres_base,nseq)) !(3,maxres_base,maxseq)
227 allocate(nres_base(3,nseq)) !(3,maxseq)
228 allocate(str_nam(nseq)) !(maxseq)
230 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),&
231 nres_base(2,i),nres_base(3,i)
232 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,&
234 ! write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
235 ! & nres_base(2,i),nres_base(3,i)
236 ! write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
240 if (weidis.eq.0.0D0) weidis=0.1D0
249 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
250 write (iout,'(a,i5)') 'nexcl: ',nexcl
251 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
253 end subroutine read_threadbase
254 !-----------------------------------------------------------------------------
256 !el subroutine read_angles(kanal,iscor,energ,iprot,*)
257 subroutine read_angles(kanal,*)
261 ! implicit real*8 (a-h,o-z)
262 ! include 'DIMENSIONS'
263 ! include 'DIMENSIONS.ZSCOPT'
264 ! include 'COMMON.INTERACT'
265 ! include 'COMMON.SBRIDGE'
266 ! include 'COMMON.GEO'
267 ! include 'COMMON.VAR'
268 ! include 'COMMON.CHAIN'
269 ! include 'COMMON.IOUNITS'
270 character(len=80) :: lineh
271 integer :: iscor,iprot,ic
272 real(kind=8) :: energ
274 subroutine read_angles(kanal,*)
279 ! implicit real*8 (a-h,o-z)
280 ! include 'DIMENSIONS'
281 ! include 'COMMON.GEO'
282 ! include 'COMMON.VAR'
283 ! include 'COMMON.CHAIN'
284 ! include 'COMMON.IOUNITS'
285 ! include 'COMMON.CONTROL'
287 ! Read angles from input
292 read(kanal,'(a80)',end=10,err=10) lineh
293 read(lineh(:5),*,err=8) ic
294 read(lineh(6:),*,err=8) energ
297 print *,'error, assuming e=1d10',lineh
301 read(lineh(18:),*,end=10,err=10) nss
303 read (lineh(20:),*,end=10,err=10) &
304 (IHPB(I),JHPB(I),I=1,NSS),iscor
306 read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
307 read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),&
310 ! print *,"energy",energ," iscor",iscor
312 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
313 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
314 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
315 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
318 ! 9/7/01 avoid 180 deg valence angle
319 if (theta(i).gt.179.99d0) theta(i)=179.99d0
321 theta(i)=deg2rad*theta(i)
322 phi(i)=deg2rad*phi(i)
323 alph(i)=deg2rad*alph(i)
324 omeg(i)=deg2rad*omeg(i)
328 end subroutine read_angles
329 !-----------------------------------------------------------------------------
330 subroutine reada(rekord,lancuch,wartosc,default)
333 character*(*) :: rekord,lancuch
334 real(kind=8) :: wartosc,default
335 integer :: iread !,ilen
337 iread=index(rekord,lancuch)
342 iread=iread+ilen(lancuch)+1
343 read (rekord(iread:),*,err=10,end=10) wartosc
348 !-----------------------------------------------------------------------------
349 subroutine readi(rekord,lancuch,wartosc,default)
352 character*(*) :: rekord,lancuch
353 integer :: wartosc,default
354 integer :: iread !,ilen
356 iread=index(rekord,lancuch)
361 iread=iread+ilen(lancuch)+1
362 read (rekord(iread:),*,err=10,end=10) wartosc
367 !-----------------------------------------------------------------------------
368 subroutine multreadi(rekord,lancuch,tablica,dim,default)
372 integer :: tablica(dim),default
373 character*(*) :: rekord,lancuch
374 character(len=80) :: aux
375 integer :: iread !,ilen
380 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
381 if (iread.eq.0) return
382 iread=iread+ilen(lancuch)+1
383 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
385 end subroutine multreadi
386 !-----------------------------------------------------------------------------
387 subroutine multreada(rekord,lancuch,tablica,dim,default)
391 real(kind=8) :: tablica(dim),default
392 character*(*) :: rekord,lancuch
393 character(len=80) :: aux
394 integer :: iread !,ilen
399 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
400 if (iread.eq.0) return
401 iread=iread+ilen(lancuch)+1
402 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
404 end subroutine multreada
405 !-----------------------------------------------------------------------------
406 subroutine card_concat(card,to_upper)
408 ! dla UNRESA to_upper jest zawsze .true.
409 ! implicit real*8 (a-h,o-z)
410 ! include 'DIMENSIONS'
411 ! include 'COMMON.IOUNITS'
413 character(len=80) :: karta !,ucase
416 read (inp,'(a)') karta
417 if (to_upper) karta=ucase(karta)
419 do while (karta(80:80).eq.'&')
420 card=card(:ilen(card)+1)//karta(:79)
421 read (inp,'(a)') karta
422 if (to_upper) karta=ucase(karta)
424 card=card(:ilen(card)+1)//karta
426 end subroutine card_concat
427 !----------------------------------------------------------------------------
428 subroutine read_afminp
431 use control_data, only:out1file
433 character*320 afmcard
436 call card_concat(afmcard,.true.)
437 call readi(afmcard,"BEG",afmbeg,0)
438 call readi(afmcard,"END",afmend,0)
439 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
440 call reada(afmcard,"VEL",velAFMconst,0.0d0)
441 print *,'FORCE=' ,forceAFMconst
442 !------ NOW PROPERTIES FOR AFM
445 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
447 distafminit=dsqrt(distafminit)
448 print *,'initdist',distafminit
450 end subroutine read_afminp
451 !-----------------------------------------------------------------------------
452 subroutine read_dist_constr
455 use geometry, only: dist
459 ! implicit real*8 (a-h,o-z)
460 ! include 'DIMENSIONS'
464 ! include 'COMMON.SETUP'
465 ! include 'COMMON.CONTROL'
466 ! include 'COMMON.CHAIN'
467 ! include 'COMMON.IOUNITS'
468 ! include 'COMMON.SBRIDGE'
469 integer,dimension(2,100) :: ifrag_,ipair_
470 real(kind=8),dimension(100) :: wfrag_,wpair_
471 character(len=640) :: controlcard
474 integer :: i,k,j,ddjk,ii,jj,itemp
475 integer :: nfrag_,npair_,ndist_
476 real(kind=8) :: dist_cut
478 ! write (iout,*) "Calling read_dist_constr"
479 ! write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
481 if(.not.allocated(ihpb)) allocate(ihpb(maxdim))
482 if(.not.allocated(jhpb)) allocate(jhpb(maxdim))
483 if(.not.allocated(dhpb)) allocate(dhpb(maxdim))
484 if(.not.allocated(dhpb1)) allocate(dhpb1(maxdim))
485 if(.not.allocated(forcon)) allocate(forcon(maxdim))
486 if(.not.allocated(fordepth)) allocate(fordepth(maxdim))
487 if ((genconstr.gt.0).and.(constr_dist.eq.11)) then
488 call gen_dist_constr2
491 call card_concat(controlcard,.true.)
492 call readi(controlcard,"NFRAG",nfrag_,0)
493 call readi(controlcard,"NPAIR",npair_,0)
494 call readi(controlcard,"NDIST",ndist_,0)
495 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
496 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
497 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
498 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
499 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
500 ! write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
501 ! write (iout,*) "IFRAG"
503 ! write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
505 ! write (iout,*) "IPAIR"
507 ! write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
509 ! if(.not.allocated(ihpb)) allocate(ihpb(maxdim))
510 ! if(.not.allocated(jhpb)) allocate(jhpb(maxdim))
511 ! if(.not.allocated(dhpb)) allocate(dhpb(maxdim))
512 ! if(.not.allocated(forcon)) allocate(forcon(maxdim))
516 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
517 if (ifrag_(2,i).gt.nstart_sup+nsup-1) &
518 ifrag_(2,i)=nstart_sup+nsup-1
519 ! write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
521 if (wfrag_(i).gt.0.0d0) then
522 do j=ifrag_(1,i),ifrag_(2,i)-1
524 ! write (iout,*) "j",j," k",k
526 if (constr_dist.eq.1) then
531 forcon(nhpb)=wfrag_(i)
532 else if (constr_dist.eq.2) then
533 if (ddjk.le.dist_cut) then
538 forcon(nhpb)=wfrag_(i)
545 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
548 if (.not.out1file .or. me.eq.king) &
549 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",&
550 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
552 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",&
553 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
560 if (wpair_(i).gt.0.0d0) then
568 do j=ifrag_(1,ii),ifrag_(2,ii)
569 do k=ifrag_(1,jj),ifrag_(2,jj)
573 forcon(nhpb)=wpair_(i)
576 if (.not.out1file .or. me.eq.king) &
577 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
578 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
580 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
581 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
588 ! read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
589 ! if (forcon(nhpb+1).gt.0.0d0) then
591 ! dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
592 if (constr_dist.eq.11) then
593 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), &
594 ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
595 fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
598 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), &
599 ibecarb(i),forcon(nhpb+1)
601 if (forcon(nhpb+1).gt.0.0d0) then
603 if (ibecarb(i).gt.0) then
607 if (dhpb(nhpb).eq.0.0d0) &
608 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
612 if (.not.out1file .or. me.eq.king) &
613 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
614 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
616 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
617 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
623 end subroutine read_dist_constr
624 !-----------------------------------------------------------------------------
625 subroutine gen_dist_constr2
628 use geometry, only: dist
633 real(kind=8) :: distance
634 if (constr_dist.eq.11) then
635 do i=nstart_sup,nstart_sup+nsup-1
636 do j=i+2,nstart_sup+nsup-1
638 if (distance.le.15.0) then
640 ihpb(nhpb)=i+nstart_seq-nstart_sup
641 jhpb(nhpb)=j+nstart_seq-nstart_sup
642 forcon(nhpb)=sqrt(0.04*distance)
643 fordepth(nhpb)=sqrt(40.0/distance)
644 dhpb(nhpb)=distance-0.1d0
645 dhpb1(nhpb)=distance+0.1d0
648 if (.not.out1file .or. me.eq.king) &
649 write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ", &
650 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
652 write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ", &
653 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
659 do i=nstart_sup,nstart_sup+nsup-1
660 do j=i+2,nstart_sup+nsup-1
662 ihpb(nhpb)=i+nstart_seq-nstart_sup
663 jhpb(nhpb)=j+nstart_seq-nstart_sup
670 end subroutine gen_dist_constr2
672 !-----------------------------------------------------------------------------
684 !-----------------------------------------------------------------------------
685 subroutine copy_to_tmp(source)
687 ! include "DIMENSIONS"
688 ! include "COMMON.IOUNITS"
689 character*(*) :: source
690 character(len=256) :: tmpfile
694 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
695 inquire(file=tmpfile,exist=ex)
697 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),&
698 " to temporary directory..."
699 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
700 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
703 end subroutine copy_to_tmp
704 !-----------------------------------------------------------------------------
705 subroutine move_from_tmp(source)
707 ! include "DIMENSIONS"
708 ! include "COMMON.IOUNITS"
709 character*(*) :: source
712 write (*,*) "Moving ",source(:ilen(source)),&
713 " from temporary directory to working directory"
714 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
715 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
717 end subroutine move_from_tmp
718 !-----------------------------------------------------------------------------
720 !-----------------------------------------------------------------------------
721 ! $Date: 1994/10/12 17:24:21 $
724 logical function find_arg(ipos,line,errflag)
726 integer, parameter :: maxlen = 80
727 character(len=80) :: line
728 character(len=1) :: empty=' ',equal='='
731 ! This function returns .TRUE., if an argument follows keyword keywd; if so
732 ! IPOS will point to the first non-blank character of the argument. Returns
733 ! .FALSE., if no argument follows the keyword; in this case IPOS points
734 ! to the first non-blank character of the next keyword.
736 do while (line(ipos:ipos) .eq. empty .and. ipos.le.maxlen)
740 if (line(ipos:ipos).eq.equal) then
743 do while (line(ipos:ipos) .eq. empty .and. ipos.le.maxlen)
746 if (ipos.gt.maxlen) errflag=.true.
752 end function find_arg
753 !-----------------------------------------------------------------------------
754 logical function find_group(iunit,jout,key1)
756 character*(*) :: key1
757 character(len=80) :: karta !,ucase
758 integer :: iunit,jout
765 do while (index(ucase(karta),key1(1:ll)).eq.0.or.lcom(1,karta))
766 read (iunit,'(a)',end=10) karta
768 write (jout,'(2a)') '> ',karta(1:78)
771 10 find_group=.false.
773 end function find_group
774 !-----------------------------------------------------------------------------
775 logical function iblnk(charc)
776 character(len=1) :: charc
779 iblnk = (n.eq.9) .or. (n.eq.10) .or. (charc.eq.' ')
782 !-----------------------------------------------------------------------------
783 integer function ilen(string)
784 character*(*) :: string
788 1 if ( ilen .gt. 0 ) then
789 if ( iblnk( string(ilen:ilen) ) ) then
796 !-----------------------------------------------------------------------------
797 integer function in_keywd_set(nkey,ikey,narg,keywd,keywdset)
798 integer :: nkey,i,ikey,narg
799 character(len=16) :: keywd,keywdset(1:nkey,0:nkey)
800 ! character(len=16) :: ucase
803 if (ucase(keywd).eq.keywdset(i,ikey)) then
809 ! No match to the allowed set of keywords if this point is reached.
812 end function in_keywd_set
813 !-----------------------------------------------------------------------------
814 character function lcase(string)
815 integer :: i, k, idiff
816 character*(*) :: string
817 character(len=1) :: c
818 character(len=40) :: chtmp
824 if (string(k+1:) .ne. ' ') then
828 idiff = ichar('a') - ichar('A')
832 if (lge(c,'A') .and. lle(c,'Z')) then
833 lcase(i:i) = char(ichar(c) + idiff)
838 !-----------------------------------------------------------------------------
839 logical function lcom(ipos,karta)
840 character(len=80) :: karta
841 character :: koment(2) = (/'!','#'/)
846 if (karta(ipos:ipos).eq.koment(i)) lcom=.true.
850 !-----------------------------------------------------------------------------
851 logical function lower_case(ch)
853 lower_case=(ch.ge.'a' .and. ch.le.'z')
855 end function lower_case
856 !-----------------------------------------------------------------------------
857 subroutine mykey(line,keywd,ipos,blankline,errflag)
859 ! This subroutine seeks a non-empty substring keywd in the string LINE.
860 ! The substring begins with the first character different from blank and
861 ! "=" encountered right to the pointer IPOS (inclusively) and terminates
862 ! at the character left to the first blank or "=". When the subroutine is
863 ! exited, the pointer IPOS is moved to the position of the terminator in LINE.
864 ! The logical variable BLANKLINE is set at .TRUE., if LINE(IPOS:) contains
865 ! only separators or the maximum length of the data line (80) has been reached.
866 ! The logical variable ERRFLAG is set at .TRUE. if the string
867 ! consists only from a "=".
868 integer, parameter :: maxlen=80
869 character(len=1) :: empty=' ',equal='=',comma=','
870 character*(*) :: keywd
871 character(len=80) :: line
872 logical :: blankline,errflag !EL,lcom
873 integer :: ipos,istart,iend
875 do while (line(ipos:ipos).eq.empty .and. (ipos.le.maxlen))
878 if (ipos.gt.maxlen .or. lcom(ipos,line) ) then
879 ! At this point the rest of the input line turned out to contain only blanks
880 ! or to be commented out.
886 ! Checks whether the current char is a separator.
887 do while (line(ipos:ipos).ne.empty .and. line(ipos:ipos).ne.equal &
888 .and. line(ipos:ipos).ne.comma .and. ipos.le.maxlen)
892 ! Error flag set to .true., if the length of the keyword was found less than 1.
893 if (iend.lt.istart) then
897 keywd=line(istart:iend)
900 !-----------------------------------------------------------------------------
901 subroutine numstr(inum,numm)
902 character(len=10) :: huj='0123456789'
903 character*(*) :: numm
904 integer :: inumm,inum,inum1,inum2
909 numm(3:3)=huj(inum2+1:inum2+1)
913 numm(2:2)=huj(inum2+1:inum2+1)
917 numm(1:1)=huj(inum2+1:inum2+1)
919 end subroutine numstr
920 !-----------------------------------------------------------------------------
921 function ucase(string)
922 integer :: i, k, idiff
923 character(*) :: string
924 character(len=len(string)) :: ucase
925 character(len=1) :: c
926 character(len=40) :: chtmp
932 if (string(k+1:) .ne. ' ') then
936 idiff = ichar('a') - ichar('A')
940 if (lge(c,'a') .and. lle(c,'z')) then
941 ucase(i:i) = char(ichar(c) - idiff)
946 !-----------------------------------------------------------------------------
948 !-----------------------------------------------------------------------------
949 subroutine pdbout(etot,tytul,iunit)
951 use geometry_data, only: c,nres
956 ! implicit real*8 (a-h,o-z)
957 ! include 'DIMENSIONS'
958 ! include 'COMMON.CHAIN'
959 ! include 'COMMON.INTERACT'
960 ! include 'COMMON.NAMES'
961 ! include 'COMMON.IOUNITS'
962 ! include 'COMMON.HEADER'
963 ! include 'COMMON.SBRIDGE'
964 ! include 'COMMON.DISTFIT'
965 ! include 'COMMON.MD'
966 !el character(len=50) :: tytul
967 character*(*) :: tytul
968 character(len=1),dimension(10) :: chainid= (/'A','B','C','D','E','F','G','H','I','J'/)
969 integer,dimension(nres) :: ica !(maxres)
972 integer :: j,iti,itj,itk,itl,i,iatom,ichain,ires,iunit
977 if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
979 write (iunit,'(3a,1pe15.5)') 'REMARK ',tytul,' ENERGY ',etot
980 !model write (iunit,'(a5,i6)') 'MODEL',1
981 if (nhfrag.gt.0) then
983 iti=itype(hfrag(1,j),1)
984 itj=itype(hfrag(2,j),1)
986 write (iunit,'(a5,i5,1x,a1,i1,2x,a3,i7,2x,a3,i7,i3,t76,i5)') &
988 restyp(iti,1),hfrag(1,j)-1,&
989 restyp(itj,1),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j)
991 write (iunit,'(a5,i5,1x,a1,i2,1x,a3,i7,2x,a3,i7,i3)') &
993 restyp(iti,1),hfrag(1,j)-1,&
994 restyp(itj,1),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j)
999 if (nbfrag.gt.0) then
1003 iti=itype(bfrag(1,j),1)
1004 itj=itype(bfrag(2,j)-1,1)
1006 write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3)') &
1008 restyp(iti,1),bfrag(1,j)-1,&
1009 restyp(itj,1),bfrag(2,j)-2,0
1011 if (bfrag(3,j).gt.bfrag(4,j)) then
1013 itk=itype(bfrag(3,j),1)
1014 itl=itype(bfrag(4,j)+1,1)
1016 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)') &
1018 restyp(itl,1),bfrag(4,j),&
1019 restyp(itk,1),bfrag(3,j)-1,-1,&
1020 "N",restyp(itk,1),bfrag(3,j)-1,&
1021 "O",restyp(iti,1),bfrag(1,j)-1
1025 itk=itype(bfrag(3,j),1)
1026 itl=itype(bfrag(4,j)-1,1)
1029 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)') &
1031 restyp(itk,1),bfrag(3,j)-1,&
1032 restyp(itl,1),bfrag(4,j)-2,1,&
1033 "N",restyp(itk,1),bfrag(3,j)-1,&
1034 "O",restyp(iti,1),bfrag(1,j)-1
1046 write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') &
1047 'SSBOND',i,'CYS',idssb(i)-nnt+1,&
1048 'CYS',jdssb(i)-nnt+1
1050 write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') &
1051 'SSBOND',i,'CYS',ihpb(i)-nnt+1-nres,&
1052 'CYS',jhpb(i)-nnt+1-nres
1063 if ((iti.eq.ntyp1).and.(iti1.eq.ntyp1)) cycle
1064 if (iti.eq.ntyp1) then
1067 write (iunit,'(a)') 'TER'
1072 write (iunit,10) iatom,restyp(iti,1),chainid(ichain),&
1073 ires,(c(j,i),j=1,3),vtot(i)
1076 write (iunit,20) iatom,restyp(iti,1),chainid(ichain),&
1077 ires,(c(j,nres+i),j=1,3),&
1082 write (iunit,'(a)') 'TER'
1084 if (itype(i,1).eq.ntyp1) cycle
1085 if (itype(i,1).eq.10 .and. itype(i+1,1).ne.ntyp1) then
1086 write (iunit,30) ica(i),ica(i+1)
1087 else if (itype(i,1).ne.10 .and. itype(i+1,1).ne.ntyp1) then
1088 write (iunit,30) ica(i),ica(i+1),ica(i)+1
1089 else if (itype(i,1).ne.10 .and. itype(i+1,1).eq.ntyp1) then
1090 write (iunit,30) ica(i),ica(i)+1
1093 if (itype(nct,1).ne.10) then
1094 write (iunit,30) ica(nct),ica(nct)+1
1098 write (iunit,30) ica(idssb(i))+1,ica(jdssb(i))+1
1100 write (iunit,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
1103 write (iunit,'(a6)') 'ENDMDL'
1104 10 FORMAT ('ATOM',I7,' CA ',A3,1X,A1,I4,4X,3F8.3,f15.3)
1105 20 FORMAT ('ATOM',I7,' CB ',A3,1X,A1,I4,4X,3F8.3,f15.3)
1106 30 FORMAT ('CONECT',8I5)
1108 end subroutine pdbout
1109 !-----------------------------------------------------------------------------
1110 subroutine MOL2out(etot,tytul)
1111 ! Prints the Cartesian coordinates of the alpha-carbons in the Tripos mol2
1113 use geometry_data, only: c
1116 ! implicit real*8 (a-h,o-z)
1117 ! include 'DIMENSIONS'
1118 ! include 'COMMON.CHAIN'
1119 ! include 'COMMON.INTERACT'
1120 ! include 'COMMON.NAMES'
1121 ! include 'COMMON.IOUNITS'
1122 ! include 'COMMON.HEADER'
1123 ! include 'COMMON.SBRIDGE'
1124 character(len=32) :: tytul,fd
1125 character(len=3) :: zahl
1126 character(len=6) :: res_num,pom !,ucase
1130 real(kind=8) :: etot
1134 #elif (defined CRAY)
1139 write (imol2,'(a)') '#'
1140 write (imol2,'(a)') &
1141 '# Creating user name: unres'
1142 write (imol2,'(2a)') '# Creation time: ',&
1144 write (imol2,'(/a)') '\@<TRIPOS>MOLECULE'
1145 write (imol2,'(a)') tytul
1146 write (imol2,'(5i5)') nct-nnt+1,nct-nnt+nss+1,nct-nnt+nss+1,0,0
1147 write (imol2,'(a)') 'SMALL'
1148 write (imol2,'(a)') 'USER_CHARGES'
1149 write (imol2,'(a)') '\@<TRIPOS>ATOM'
1151 write (zahl,'(i3)') i
1152 pom=ucase(restyp(itype(i,1),1))
1153 res_num = pom(:3)//zahl(2:)
1154 write (imol2,10) i,(c(j,i),j=1,3),i,res_num,0.0
1156 write (imol2,'(a)') '\@<TRIPOS>BOND'
1158 write (imol2,'(i5,2i6,i2)') i-nnt+1,i-nnt+1,i-nnt+2,1
1161 write (imol2,'(i5,2i6,i2)') nct-nnt+i,ihpb(i),jhpb(i),1
1163 write (imol2,'(a)') '\@<TRIPOS>SUBSTRUCTURE'
1165 write (zahl,'(i3)') i
1166 pom = ucase(restyp(itype(i,1),1))
1167 res_num = pom(:3)//zahl(2:)
1168 write (imol2,30) i-nnt+1,res_num,i-nnt+1,0
1170 10 FORMAT (I7,' CA ',3F10.4,' C.3',I8,1X,A,F11.4,' ****')
1171 30 FORMAT (I7,1x,A,I14,' RESIDUE',I13,' **** ****')
1173 end subroutine MOL2out
1174 !-----------------------------------------------------------------------------
1178 use energy_data, only: itype
1180 ! implicit real*8 (a-h,o-z)
1181 ! include 'DIMENSIONS'
1182 ! include 'COMMON.IOUNITS'
1183 ! include 'COMMON.CHAIN'
1184 ! include 'COMMON.VAR'
1185 ! include 'COMMON.LOCAL'
1186 ! include 'COMMON.INTERACT'
1187 ! include 'COMMON.NAMES'
1188 ! include 'COMMON.GEO'
1189 ! include 'COMMON.TORSION'
1193 write (iout,'(/a)') 'Geometry of the virtual chain.'
1194 write (iout,'(7a)') ' Res ',' d',' Theta',&
1195 ' Phi',' Dsc',' Alpha',' Omega'
1198 write (iout,'(a3,i4,6f10.3)') restyp(iti,1),i,vbld(i),&
1199 rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),rad2deg*alph(i),&
1203 end subroutine intout
1204 !-----------------------------------------------------------------------------
1206 subroutine briefout(it,ener,free)!,plik)
1208 subroutine briefout(it,ener)
1213 ! implicit real*8 (a-h,o-z)
1214 ! include 'DIMENSIONS'
1215 ! include 'COMMON.IOUNITS'
1216 ! include 'COMMON.CHAIN'
1217 ! include 'COMMON.VAR'
1218 ! include 'COMMON.LOCAL'
1219 ! include 'COMMON.INTERACT'
1220 ! include 'COMMON.NAMES'
1221 ! include 'COMMON.GEO'
1222 ! include 'COMMON.SBRIDGE'
1223 ! print '(a,i5)',intname,igeom
1226 real(kind=8) :: ener,free
1227 ! character(len=80) :: plik
1230 #if defined(AIX) || defined(PGI)
1231 open (igeom,file=intname,position='append')
1233 open (igeom,file=intname,access='append')
1242 WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,NSS)
1244 WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,8)
1246 WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,NSS)
1248 WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,9)
1250 WRITE (igeom,190) (IHPB(I),JHPB(I),I=10,NSS)
1252 ! IF (nvar.gt.nphi) WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
1253 WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
1254 WRITE (igeom,200) (RAD2DEG*PHI(I),I=4,NRES)
1255 ! if (nvar.gt.nphi+ntheta) then
1256 write (igeom,200) (rad2deg*alph(i),i=2,nres-1)
1257 write (igeom,200) (rad2deg*omeg(i),i=2,nres-1)
1260 180 format (I5,F12.3,I2,9(1X,2I3))
1261 190 format (3X,11(1X,2I3))
1264 end subroutine briefout
1265 !-----------------------------------------------------------------------------
1267 subroutine fdate(fd)
1268 character(len=32) :: fd
1271 end subroutine fdate
1273 !-----------------------------------------------------------------------------
1275 real(kind=8) function gyrate(jcon)
1277 real(kind=8) function gyrate()
1280 use geometry_data, only: c
1283 ! implicit real*8 (a-h,o-z)
1284 ! include 'DIMENSIONS'
1285 ! include 'COMMON.INTERACT'
1286 ! include 'COMMON.CHAIN'
1288 real(kind=8),dimension(3) :: cen
1298 cen(j)=cen(j)+c(j,i)
1302 cen(j)=cen(j)/dble(nct-nnt+1)
1307 rg = rg + (c(j,i)-cen(j))**2
1311 gyrate = dsqrt(rg/dble(nct-nnt+1))
1313 gyrate = sqrt(rg/dble(nct-nnt+1))
1318 !-----------------------------------------------------------------------------
1320 subroutine reads(rekord,lancuch,wartosc,default)
1322 character*(*) :: rekord,lancuch,wartosc,default
1323 character(len=80) :: aux
1324 integer :: lenlan,lenrec,iread,ireade
1328 lenlan=ilen(lancuch)
1330 iread=index(rekord,lancuch(:lenlan)//"=")
1331 ! print *,"rekord",rekord," lancuch",lancuch
1332 ! print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
1333 if (iread.eq.0) then
1337 iread=iread+lenlan+1
1338 ! print *,"iread",iread
1339 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
1340 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
1342 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
1344 ! print *,"iread",iread
1345 if (iread.gt.lenrec) then
1350 ! print *,"ireade",ireade
1351 do while (ireade.lt.lenrec .and. &
1352 .not.iblnk(rekord(ireade:ireade)))
1355 wartosc=rekord(iread:ireade)
1357 end subroutine reads
1359 !-----------------------------------------------------------------------------
1361 !-----------------------------------------------------------------------------
1362 subroutine permut(isym)
1364 use geometry_data, only: tabperm
1365 ! implicit real*8 (a-h,o-z)
1366 ! include 'DIMENSIONS'
1367 ! include 'COMMON.LOCAL'
1368 ! include 'COMMON.VAR'
1369 ! include 'COMMON.CHAIN'
1370 ! include 'COMMON.INTERACT'
1371 ! include 'COMMON.IOUNITS'
1372 ! include 'COMMON.GEO'
1373 ! include 'COMMON.NAMES'
1374 ! include 'COMMON.CONTROL'
1379 integer,dimension(isym) :: a
1380 ! parameter(n=symetr)
1393 10 print *,(a(i),i=1,n)
1397 ! write (iout,*) "tututu", kkk
1399 if(nextp(n,a)) go to 10
1401 end subroutine permut
1402 !-----------------------------------------------------------------------------
1403 logical function nextp(n,a)
1405 integer :: n,i,j,k,t
1407 integer,dimension(n) :: a
1409 10 if(a(i).lt.a(i+1)) go to 20
1426 if(a(j).lt.a(i)) go to 40
1433 !-----------------------------------------------------------------------------
1434 !-----------------------------------------------------------------------------