2 !-----------------------------------------------------------------------
6 !-----------------------------------------------------------------------------
7 ! Max. number of AA residues
8 integer,parameter :: maxres=101000!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 !-----------------------------------------------------------------------------
29 subroutine read_bridge
30 ! Read information about disulfide bridges.
31 use geometry_data, only: nres
33 use control_data, only:out1file
35 ! implicit real*8 (a-h,o-z)
36 ! include 'DIMENSIONS'
40 ! include 'COMMON.IOUNITS'
41 ! include 'COMMON.GEO'
42 ! include 'COMMON.VAR'
43 ! include 'COMMON.INTERACT'
44 ! include 'COMMON.LOCAL'
45 ! include 'COMMON.NAMES'
46 ! include 'COMMON.CHAIN'
47 ! include 'COMMON.FFIELD'
48 ! include 'COMMON.SBRIDGE'
49 ! include 'COMMON.HEADER'
50 ! include 'COMMON.CONTROL'
51 ! include 'COMMON.DBASE'
52 ! include 'COMMON.THREAD'
53 ! include 'COMMON.TIME1'
54 ! include 'COMMON.SETUP'
58 ! Read bridging residues.
64 read (inp,*) (iss(i),i=1,ns)
68 if(me.eq.king.or..not.out1file) &
69 write (iout,*) 'ns=',ns
71 write(iout,*) ' iss:',(iss(i),i=1,ns)
72 ! Check whether the specified bridging residues are cystines.
74 write(iout,*) i,iss(i)
75 if (itype(iss(i),1).ne.1) then
76 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)') &
77 'Do you REALLY think that the residue ',&
78 restyp(itype(iss(i),1),1),i,&
79 ' can form a disulfide bridge?!!!'
80 write (*,'(2a,i3,a)') &
81 'Do you REALLY think that the residue ',&
82 restyp(itype(iss(i),1),1),i,&
83 ' can form a disulfide bridge?!!!'
85 call MPI_Finalize(MPI_COMM_WORLD,ierror)
91 if(.not.allocated(ihpb)) allocate(ihpb(ns))
92 if(.not.allocated(jhpb)) allocate(jhpb(ns))
94 ! Read preformed bridges.
98 if(.not.allocated(ihpb)) allocate(ihpb(nss))
99 if(.not.allocated(jhpb)) allocate(jhpb(nss))
100 read (inp,*) (ihpb(i),jhpb(i),i=1,nss)
102 if(.not.allocated(dhpb)) allocate(dhpb(nss))
103 if(.not.allocated(forcon)) allocate(forcon(nss))!(maxdim) !el maxdim=(maxres-1)*(maxres-2)/2
104 if(.not.allocated(dhpb1)) allocate(dhpb1(nss))
105 if(.not.allocated(ibecarb)) allocate(ibecarb(nss))
106 ! el Initialize the bridge array
111 !--------------------
113 write(iout,*)'nss=',nss !el,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
115 write(iout,*)'ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
117 ! Check if the residues involved in bridges are in the specified list of
121 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j) &
122 .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
123 write (iout,'(a,i3,a)') 'Disulfide pair',i,&
124 ' contains residues present in other pairs.'
125 write (*,'(a,i3,a)') 'Disulfide pair',i,&
126 ' contains residues present in other pairs.'
128 call MPI_Finalize(MPI_COMM_WORLD,ierror)
134 if (ihpb(i).eq.iss(j)) goto 10
136 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
139 if (jhpb(i).eq.iss(j)) goto 20
141 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
152 write(iout,*) "end read_bridge"
154 end subroutine read_bridge
155 !-----------------------------------------------------------------------------
156 subroutine read_x(kanal,*)
160 use geometry, only:int_from_cart1
161 ! implicit real*8 (a-h,o-z)
162 ! include 'DIMENSIONS'
163 ! include 'COMMON.GEO'
164 ! include 'COMMON.VAR'
165 ! include 'COMMON.CHAIN'
166 ! include 'COMMON.IOUNITS'
167 ! include 'COMMON.CONTROL'
168 ! include 'COMMON.LOCAL'
169 ! include 'COMMON.INTERACT'
170 ! Read coordinates from input
173 integer :: l,k,j,i,kanal
175 read(kanal,'(8f10.5)',end=10,err=10) &
176 ((c(l,k),l=1,3),k=1,nres),&
177 ((c(l,k+nres),l=1,3),k=nnt,nct)
180 c(j,2*nres)=c(j,nres)
182 call int_from_cart1(.false.)
185 dc(j,i)=c(j,i+1)-c(j,i)
186 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
190 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
192 dc(j,i+nres)=c(j,i+nres)-c(j,i)
193 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
200 end subroutine read_x
201 !-----------------------------------------------------------------------------
202 subroutine read_threadbase
207 ! implicit real*8 (a-h,o-z)
208 ! include 'DIMENSIONS'
209 ! include 'COMMON.IOUNITS'
210 ! include 'COMMON.GEO'
211 ! include 'COMMON.VAR'
212 ! include 'COMMON.INTERACT'
213 ! include 'COMMON.LOCAL'
214 ! include 'COMMON.NAMES'
215 ! include 'COMMON.CHAIN'
216 ! include 'COMMON.FFIELD'
217 ! include 'COMMON.SBRIDGE'
218 ! include 'COMMON.HEADER'
219 ! include 'COMMON.CONTROL'
220 ! include 'COMMON.DBASE'
221 ! include 'COMMON.THREAD'
222 ! include 'COMMON.TIME1'
227 ! Read pattern database for threading.
229 allocate(cart_base(3,maxres_base,nseq)) !(3,maxres_base,maxseq)
230 allocate(nres_base(3,nseq)) !(3,maxseq)
231 allocate(str_nam(nseq)) !(maxseq)
233 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),&
234 nres_base(2,i),nres_base(3,i)
235 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,&
237 ! write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
238 ! & nres_base(2,i),nres_base(3,i)
239 ! write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
243 if (weidis.eq.0.0D0) weidis=0.1D0
252 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
253 write (iout,'(a,i5)') 'nexcl: ',nexcl
254 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
256 end subroutine read_threadbase
257 !-----------------------------------------------------------------------------
259 !el subroutine read_angles(kanal,iscor,energ,iprot,*)
260 subroutine read_angles(kanal,*)
264 ! implicit real*8 (a-h,o-z)
265 ! include 'DIMENSIONS'
266 ! include 'DIMENSIONS.ZSCOPT'
267 ! include 'COMMON.INTERACT'
268 ! include 'COMMON.SBRIDGE'
269 ! include 'COMMON.GEO'
270 ! include 'COMMON.VAR'
271 ! include 'COMMON.CHAIN'
272 ! include 'COMMON.IOUNITS'
273 character(len=80) :: lineh
274 integer :: iscor,iprot,ic
275 real(kind=8) :: energ
277 subroutine read_angles(kanal,*)
283 ! implicit real*8 (a-h,o-z)
284 ! include 'DIMENSIONS'
285 ! include 'COMMON.GEO'
286 ! include 'COMMON.VAR'
287 ! include 'COMMON.CHAIN'
288 ! include 'COMMON.IOUNITS'
289 ! include 'COMMON.CONTROL'
291 ! Read angles from input
296 read(kanal,'(a80)',end=10,err=10) lineh
297 read(lineh(:5),*,err=8) ic
298 read(lineh(6:),*,err=8) energ
301 print *,'error, assuming e=1d10',lineh
305 read(lineh(18:),*,end=10,err=10) nss
307 read (lineh(20:),*,end=10,err=10) &
308 (IHPB(I),JHPB(I),I=1,NSS),iscor
310 read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
311 read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),&
314 ! print *,"energy",energ," iscor",iscor
316 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
317 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
318 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
319 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
322 ! 9/7/01 avoid 180 deg valence angle
323 if (theta(i).gt.179.99d0) theta(i)=179.99d0
325 theta(i)=deg2rad*theta(i)
326 phi(i)=deg2rad*phi(i)
327 alph(i)=deg2rad*alph(i)
328 omeg(i)=deg2rad*omeg(i)
332 end subroutine read_angles
333 !-----------------------------------------------------------------------------
334 subroutine reada(rekord,lancuch,wartosc,default)
337 character*(*) :: rekord,lancuch
338 real(kind=8) :: wartosc,default
339 integer :: iread !,ilen
341 iread=index(rekord,lancuch)
346 iread=iread+ilen(lancuch)+1
347 read (rekord(iread:),*,err=10,end=10) wartosc
352 !-----------------------------------------------------------------------------
353 subroutine readi(rekord,lancuch,wartosc,default)
356 character*(*) :: rekord,lancuch
357 integer :: wartosc,default
358 integer :: iread !,ilen
360 iread=index(rekord,lancuch)
365 iread=iread+ilen(lancuch)+1
366 read (rekord(iread:),*,err=10,end=10) wartosc
371 !-----------------------------------------------------------------------------
372 subroutine multreadi(rekord,lancuch,tablica,dim,default)
376 integer :: tablica(dim),default
377 character*(*) :: rekord,lancuch
378 character(len=80) :: aux
379 integer :: iread !,ilen
384 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
385 if (iread.eq.0) return
386 iread=iread+ilen(lancuch)+1
387 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
389 end subroutine multreadi
390 !-----------------------------------------------------------------------------
391 subroutine multreada(rekord,lancuch,tablica,dim,default)
395 real(kind=8) :: tablica(dim),default
396 character*(*) :: rekord,lancuch
397 character(len=80) :: aux
398 integer :: iread !,ilen
403 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
404 if (iread.eq.0) return
405 iread=iread+ilen(lancuch)+1
406 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
408 end subroutine multreada
409 !-----------------------------------------------------------------------------
410 subroutine card_concat(card,to_upper)
412 ! dla UNRESA to_upper jest zawsze .true.
413 ! implicit real*8 (a-h,o-z)
414 ! include 'DIMENSIONS'
415 ! include 'COMMON.IOUNITS'
417 character(len=80) :: karta !,ucase
420 read (inp,'(a)') karta
421 if (to_upper) karta=ucase(karta)
423 do while (karta(80:80).eq.'&')
424 card=card(:ilen(card)+1)//karta(:79)
425 read (inp,'(a)') karta
426 if (to_upper) karta=ucase(karta)
428 card=card(:ilen(card)+1)//karta
430 end subroutine card_concat
431 !----------------------------------------------------------------------------
432 subroutine read_afminp
435 use control_data, only:out1file
437 character*320 afmcard
438 real, dimension(3) ::cbeg,cend
441 call card_concat(afmcard,.true.)
442 call readi(afmcard,"BEG",afmbeg,0)
443 call readi(afmcard,"END",afmend,0)
444 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
445 call reada(afmcard,"VEL",velAFMconst,0.0d0)
446 print *,'FORCE=' ,forceAFMconst
449 if (afmbeg.eq.-1) then
450 read(inp,*) nbegafmmat,(afmbegcentr(i),i=1,nbegafmmat)
453 cbeg(j)=cbeg(j)+c(j,afmbegcentr(i))/nbegafmmat
461 if (afmend.eq.-1) then
462 read(inp,*) nendafmmat,(afmendcentr(i),i=1,nendafmmat)
465 cend(j)=cend(j)+c(j,afmendcentr(i))/nendafmmat
471 !------ NOW PROPERTIES FOR AFM
474 distafminit=(cend(i)-cbeg(i))**2+distafminit
476 distafminit=dsqrt(distafminit)
477 print *,'initdist',distafminit
479 end subroutine read_afminp
480 !C------------------------------------------------------------
481 subroutine read_afmnano
484 use control_data, only:out1file
488 read(inp,*,err=11) inanomove,(inanotab(i),i=1,inanomove),forcenanoconst
491 cm=cm+c(3,inanotab(i))
494 distnanoinit=cm-tubecenter(3)
498 ", number of center, their index and forceconstance"
501 end subroutine read_afmnano
503 !-----------------------------------------------------------------------------
504 subroutine read_dist_constr
507 use geometry, only: dist
511 ! implicit real*8 (a-h,o-z)
512 ! include 'DIMENSIONS'
516 ! include 'COMMON.SETUP'
517 ! include 'COMMON.CONTROL'
518 ! include 'COMMON.CHAIN'
519 ! include 'COMMON.IOUNITS'
520 ! include 'COMMON.SBRIDGE'
521 integer,dimension(2,100) :: ifrag_,ipair_
522 real(kind=8),dimension(100) :: wfrag_,wpair_
523 character(len=640) :: controlcard
526 integer :: i,k,j,ii,jj,itemp,mnumkk,mnumjj,k1,j1
527 integer :: nfrag_,npair_,ndist_
528 real(kind=8) :: dist_cut,ddjk
530 ! write (iout,*) "Calling read_dist_constr"
531 ! write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
533 if(.not.allocated(ihpb)) allocate(ihpb(maxdim))
534 if(.not.allocated(jhpb)) allocate(jhpb(maxdim))
535 if(.not.allocated(dhpb)) allocate(dhpb(maxdim))
536 if(.not.allocated(dhpb1)) allocate(dhpb1(maxdim))
537 if(.not.allocated(forcon)) allocate(forcon(maxdim))
538 if(.not.allocated(fordepth)) allocate(fordepth(maxdim))
539 if(.not.allocated(ibecarb)) allocate(ibecarb(maxdim))
540 if ((genconstr.gt.0).and.(constr_dist.eq.11)) then
541 call gen_dist_constr2
544 call card_concat(controlcard,.true.)
545 call readi(controlcard,"NFRAG",nfrag_,0)
546 call readi(controlcard,"NPAIR",npair_,0)
547 call readi(controlcard,"NDIST",ndist_,0)
548 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
549 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
550 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
551 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
552 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
553 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
554 ! write (iout,*) "IFRAG"
556 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
558 ! write (iout,*) "IPAIR"
560 ! write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
562 ! if(.not.allocated(ihpb)) allocate(ihpb(maxdim))
563 ! if(.not.allocated(jhpb)) allocate(jhpb(maxdim))
564 ! if(.not.allocated(dhpb)) allocate(dhpb(maxdim))
565 ! if(.not.allocated(forcon)) allocate(forcon(maxdim))
569 ! if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
570 ! if (ifrag_(2,i).gt.nstart_sup+nsup-1) &
571 ! ifrag_(2,i)=nstart_sup+nsup-1
572 ! write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
574 if (wfrag_(i).gt.0.0d0) then
575 do j=ifrag_(1,i),ifrag_(2,i)-1
577 write (iout,*) "j",j," k",k,nres
580 if (j.gt.nres) j1=j-nres
581 if (k.gt.nres) k1=k-nres
585 if ((itype(k1,mnumkk).gt.ntyp_molec(mnumkk)).or.&
586 (itype(j1,mnumjj).gt.ntyp_molec(mnumjj))) cycle
587 write (iout,*) "j",j," k",k,itype(k1,mnumkk),itype(j1,mnumjj)
589 if (constr_dist.eq.1) then
594 forcon(nhpb)=wfrag_(i)
595 else if (constr_dist.eq.2) then
596 if (ddjk.le.dist_cut) then
601 forcon(nhpb)=wfrag_(i)
608 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
611 if (.not.out1file .or. me.eq.king) &
612 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",&
613 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
615 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",&
616 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
623 if (wpair_(i).gt.0.0d0) then
631 do j=ifrag_(1,ii),ifrag_(2,ii)
632 do k=ifrag_(1,jj),ifrag_(2,jj)
636 forcon(nhpb)=wpair_(i)
639 if (.not.out1file .or. me.eq.king) &
640 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
641 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
643 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
644 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
651 ! read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
652 ! if (forcon(nhpb+1).gt.0.0d0) then
654 ! dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
655 if (constr_dist.eq.11) then
656 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), &
657 ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
658 !EL fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
659 fordepth(nhpb+1)=fordepth(nhpb+1)**(0.25d0)
660 forcon(nhpb+1)=forcon(nhpb+1)**(0.25d0)
663 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), &
664 ibecarb(i),forcon(nhpb+1)
666 if (forcon(nhpb+1).gt.0.0d0) then
668 if (ibecarb(i).gt.0) then
672 if (dhpb(nhpb).eq.0.0d0) &
673 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
677 if (.not.out1file .or. me.eq.king) &
678 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
679 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
681 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
682 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
688 end subroutine read_dist_constr
689 !-----------------------------------------------------------------------------
690 subroutine gen_dist_constr2
693 use geometry, only: dist
698 real(kind=8) :: distance
699 if (constr_dist.eq.11) then
700 do i=nstart_sup,nstart_sup+nsup-1
701 do j=i+2,nstart_sup+nsup-1
703 if (distance.le.15.0) then
705 ihpb(nhpb)=i+nstart_seq-nstart_sup
706 jhpb(nhpb)=j+nstart_seq-nstart_sup
707 forcon(nhpb)=sqrt(0.04*distance)
708 fordepth(nhpb)=sqrt(40.0/distance)
709 dhpb(nhpb)=distance-0.1d0
710 dhpb1(nhpb)=distance+0.1d0
713 if (.not.out1file .or. me.eq.king) &
714 write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ", &
715 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
717 write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ", &
718 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
724 do i=nstart_sup,nstart_sup+nsup-1
725 do j=i+2,nstart_sup+nsup-1
727 ihpb(nhpb)=i+nstart_seq-nstart_sup
728 jhpb(nhpb)=j+nstart_seq-nstart_sup
735 end subroutine gen_dist_constr2
737 !-----------------------------------------------------------------------------
749 !-----------------------------------------------------------------------------
750 subroutine copy_to_tmp(source)
752 ! include "DIMENSIONS"
753 ! include "COMMON.IOUNITS"
754 character*(*) :: source
755 character(len=256) :: tmpfile
759 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
760 inquire(file=tmpfile,exist=ex)
762 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),&
763 " to temporary directory..."
764 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
765 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
768 end subroutine copy_to_tmp
769 !-----------------------------------------------------------------------------
770 subroutine move_from_tmp(source)
772 ! include "DIMENSIONS"
773 ! include "COMMON.IOUNITS"
774 character*(*) :: source
777 write (*,*) "Moving ",source(:ilen(source)),&
778 " from temporary directory to working directory"
779 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
780 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
782 end subroutine move_from_tmp
783 !-----------------------------------------------------------------------------
785 !-----------------------------------------------------------------------------
786 ! $Date: 1994/10/12 17:24:21 $
789 logical function find_arg(ipos,line,errflag)
791 integer, parameter :: maxlen = 80
792 character(len=80) :: line
793 character(len=1) :: empty=' ',equal='='
796 ! This function returns .TRUE., if an argument follows keyword keywd; if so
797 ! IPOS will point to the first non-blank character of the argument. Returns
798 ! .FALSE., if no argument follows the keyword; in this case IPOS points
799 ! to the first non-blank character of the next keyword.
801 do while (line(ipos:ipos) .eq. empty .and. ipos.le.maxlen)
805 if (line(ipos:ipos).eq.equal) then
808 do while (line(ipos:ipos) .eq. empty .and. ipos.le.maxlen)
811 if (ipos.gt.maxlen) errflag=.true.
817 end function find_arg
818 !-----------------------------------------------------------------------------
819 logical function find_group(iunit,jout,key1)
821 character*(*) :: key1
822 character(len=80) :: karta !,ucase
823 integer :: iunit,jout
830 do while (index(ucase(karta),key1(1:ll)).eq.0.or.lcom(1,karta))
831 read (iunit,'(a)',end=10) karta
833 write (jout,'(2a)') '> ',karta(1:78)
836 10 find_group=.false.
838 end function find_group
840 !-----------------------------------------------------------------------------
841 logical function iblnk(charc)
842 character(len=1) :: charc
845 iblnk = (n.eq.9) .or. (n.eq.10) .or. (charc.eq.' ')
848 !-----------------------------------------------------------------------------
850 integer function ilen(string)
851 character*(*) :: string
855 1 if ( ilen .gt. 0 ) then
856 if ( iblnk( string(ilen:ilen) ) ) then
863 !-----------------------------------------------------------------------------
864 integer function in_keywd_set(nkey,ikey,narg,keywd,keywdset)
865 integer :: nkey,i,ikey,narg
866 character(len=16) :: keywd,keywdset(1:nkey,0:nkey)
867 ! character(len=16) :: ucase
870 if (ucase(keywd).eq.keywdset(i,ikey)) then
876 ! No match to the allowed set of keywords if this point is reached.
879 end function in_keywd_set
880 !-----------------------------------------------------------------------------
881 character function lcase(string)
882 integer :: i, k, idiff
883 character*(*) :: string
884 character(len=1) :: c
885 character(len=40) :: chtmp
891 if (string(k+1:) .ne. ' ') then
895 idiff = ichar('a') - ichar('A')
899 if (lge(c,'A') .and. lle(c,'Z')) then
900 lcase(i:i) = char(ichar(c) + idiff)
905 !-----------------------------------------------------------------------------
906 logical function lcom(ipos,karta)
907 character(len=80) :: karta
908 character :: koment(2) = (/'!','#'/)
913 if (karta(ipos:ipos).eq.koment(i)) lcom=.true.
917 !-----------------------------------------------------------------------------
918 logical function lower_case(ch)
920 lower_case=(ch.ge.'a' .and. ch.le.'z')
922 end function lower_case
923 !-----------------------------------------------------------------------------
924 subroutine mykey(line,keywd,ipos,blankline,errflag)
926 ! This subroutine seeks a non-empty substring keywd in the string LINE.
927 ! The substring begins with the first character different from blank and
928 ! "=" encountered right to the pointer IPOS (inclusively) and terminates
929 ! at the character left to the first blank or "=". When the subroutine is
930 ! exited, the pointer IPOS is moved to the position of the terminator in LINE.
931 ! The logical variable BLANKLINE is set at .TRUE., if LINE(IPOS:) contains
932 ! only separators or the maximum length of the data line (80) has been reached.
933 ! The logical variable ERRFLAG is set at .TRUE. if the string
934 ! consists only from a "=".
935 integer, parameter :: maxlen=80
936 character(len=1) :: empty=' ',equal='=',comma=','
937 character*(*) :: keywd
938 character(len=80) :: line
939 logical :: blankline,errflag !EL,lcom
940 integer :: ipos,istart,iend
942 do while (line(ipos:ipos).eq.empty .and. (ipos.le.maxlen))
945 if (ipos.gt.maxlen .or. lcom(ipos,line) ) then
946 ! At this point the rest of the input line turned out to contain only blanks
947 ! or to be commented out.
953 ! Checks whether the current char is a separator.
954 do while (line(ipos:ipos).ne.empty .and. line(ipos:ipos).ne.equal &
955 .and. line(ipos:ipos).ne.comma .and. ipos.le.maxlen)
959 ! Error flag set to .true., if the length of the keyword was found less than 1.
960 if (iend.lt.istart) then
964 keywd=line(istart:iend)
967 !-----------------------------------------------------------------------------
968 subroutine numstr(inum,numm)
969 character(len=10) :: huj='0123456789'
970 character*(*) :: numm
971 integer :: inumm,inum,inum1,inum2
976 numm(3:3)=huj(inum2+1:inum2+1)
980 numm(2:2)=huj(inum2+1:inum2+1)
984 numm(1:1)=huj(inum2+1:inum2+1)
986 end subroutine numstr
988 !-----------------------------------------------------------------------------
989 function ucase(string)
990 integer :: i, k, idiff
991 character(*) :: string
992 character(len=len(string)) :: ucase
993 character(len=1) :: c
994 character(len=40) :: chtmp
1000 if (string(k+1:) .ne. ' ') then
1004 idiff = ichar('a') - ichar('A')
1008 if (lge(c,'a') .and. lle(c,'z')) then
1009 ucase(i:i) = char(ichar(c) - idiff)
1014 !-----------------------------------------------------------------------------
1016 !-----------------------------------------------------------------------------
1017 subroutine pdbout(etot,tytul,iunit)
1019 use geometry_data, only: c,nres,boxxsize,boxysize,boxzsize
1021 ! use energy, only: to_box
1025 ! implicit real*8 (a-h,o-z)
1026 ! include 'DIMENSIONS'
1027 ! include 'COMMON.CHAIN'
1028 ! include 'COMMON.INTERACT'
1029 ! include 'COMMON.NAMES'
1030 ! include 'COMMON.IOUNITS'
1031 ! include 'COMMON.HEADER'
1032 ! include 'COMMON.SBRIDGE'
1033 ! include 'COMMON.DISTFIT'
1034 ! include 'COMMON.MD'
1035 !el character(len=50) :: tytul
1036 character*(*) :: tytul
1037 character(len=1),dimension(58) :: chainid= (/'A','B','C','D','E','F','G','H','I','J',&
1038 'K','L','M','O','Q','P','R','S','T','U','W','X','Y','Z','a','b','c','d','e','f',&
1039 'g','h','i','j','k','l','m','n','o','p','r','s','t','u','w','x','y','z',&
1040 '0','1','2','3','4','5','6','7','8','9'/)
1042 integer,dimension(maxres) :: ica !(maxres)
1044 integer,dimension(nres) :: ica !(maxres)
1048 integer :: j,iti,itj,itk,itl,i,iatom,ichain,ires,iunit
1049 real(kind=8) :: etot,xi,yi,zi
1053 if(.not.allocated(molnum)) allocate(molnum(maxres))
1054 ! molnum(:)=mnumi(:)
1055 if(.not.allocated(vtot)) allocate(vtot(maxres*2)) !(maxres2)
1056 if(.not.allocated(istype)) allocate(istype(maxres))
1058 if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
1061 write (iunit,'(3a,1pe15.5)') 'REMARK ',tytul,' ENERGY ',etot
1062 !model write (iunit,'(a5,i6)') 'MODEL',1
1064 if (nhfrag.gt.0) then
1066 iti=itype(hfrag(1,j),1)
1067 itj=itype(hfrag(2,j),1)
1069 write (iunit,'(a5,i5,1x,a1,i1,2x,a3,i7,2x,a3,i7,i3,t76,i5)') &
1071 restyp(iti,1),hfrag(1,j)-1,&
1072 restyp(itj,1),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j)
1074 write (iunit,'(a5,i5,1x,a1,i2,1x,a3,i7,2x,a3,i7,i3)') &
1076 restyp(iti,1),hfrag(1,j)-1,&
1077 restyp(itj,1),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j)
1082 if (nbfrag.gt.0) then
1086 iti=itype(bfrag(1,j),1)
1087 itj=itype(bfrag(2,j)-1,1)
1089 write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3)') &
1091 restyp(iti,1),bfrag(1,j)-1,&
1092 restyp(itj,1),bfrag(2,j)-2,0
1094 if (bfrag(3,j).gt.bfrag(4,j)) then
1096 itk=itype(bfrag(3,j),1)
1097 itl=itype(bfrag(4,j)+1,1)
1099 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)') &
1101 restyp(itl,1),bfrag(4,j),&
1102 restyp(itk,1),bfrag(3,j)-1,-1,&
1103 "N",restyp(itk,1),bfrag(3,j)-1,&
1104 "O",restyp(iti,1),bfrag(1,j)-1
1108 itk=itype(bfrag(3,j),1)
1109 itl=itype(bfrag(4,j)-1,1)
1112 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)') &
1114 restyp(itk,1),bfrag(3,j)-1,&
1115 restyp(itl,1),bfrag(4,j)-2,1,&
1116 "N",restyp(itk,1),bfrag(3,j)-1,&
1117 "O",restyp(iti,1),bfrag(1,j)-1
1129 write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') &
1130 'SSBOND',i,'CYS',idssb(i)-nnt+1,&
1131 'CYS',jdssb(i)-nnt+1
1133 write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') &
1134 'SSBOND',i,'CYS',ihpb(i)-nnt+1-nres,&
1135 'CYS',jhpb(i)-nnt+1-nres
1143 write(iout,*) "TUTUT"
1145 write(iout,*), "coord",c(1,i),c(2,i),c(3,i)
1146 iti=itype(i,molnum(i))
1148 if (molnum(i+1).eq.0) then
1149 iti1=ntyp1_molec(molnum(i))
1151 iti1=itype(i+1,molnum(i+1))
1153 if ((iti.eq.ntyp1_molec(molnum(i))).and.(iti1.eq.ntyp1_molec(molnum(i)))) cycle
1155 if (iti.eq.ntyp1_molec(molnum(i)).and.(molnum(i+1).eq.5)) cycle
1157 if (iti.eq.ntyp1_molec(molnum(i))) then
1159 if (ichain.gt.58) ichain=1
1161 write (iunit,'(a)') 'TER'
1166 if (molnum(i).eq.1) then
1168 write (iunit,10) iatom,restyp(iti,molnum(i)),chainid(ichain),&
1169 ires,(c(j,i),j=1,3),vtot(i)
1170 elseif(molnum(i).eq.2) then
1171 if (istype(i).eq.0) istype(i)=1
1172 write (iunit,40) iatom,sugartyp(istype(i)),restyp(iti,2), &
1173 chainid(ichain),ires,(c(j,i),j=1,3),vtot(i)
1174 elseif (molnum(i).eq.4) then
1178 ! xi=dmod(xi,boxxsize)
1179 ! if (xi.lt.0.0d0) xi=xi+boxxsize
1180 ! yi=dmod(yi,boxysize)
1181 ! if (yi.lt.0.0d0) yi=yi+boxysize
1182 ! zi=dmod(zi,boxzsize)
1183 ! if (zi.lt.0.0d0) zi=zi+boxzsize
1184 write (iunit,60) iatom,restyp(iti,molnum(i)),chainid(ichain),&
1185 ires,xi,yi,zi,vtot(i)
1187 write (iunit,60) iatom,restyp(iti,molnum(i)),chainid(ichain),&
1188 ires,(c(j,i),j=1,3),vtot(i)
1190 if ((iti.ne.10).and.(molnum(i).ne.5)) then
1192 if (molnum(i).eq.1) then
1193 write (iunit,20) iatom,restyp(iti,1),chainid(ichain),&
1194 ires,(c(j,nres+i),j=1,3),&
1196 else if (molnum(i).eq.2) then
1197 if (istype(i).eq.0) istype(i)=1
1198 write (iunit,50) iatom,sugartyp(istype(i)),restyp(iti,2), &
1199 chainid(ichain),ires,(c(j,nres+i),j=1,3),vtot(i+nres)
1205 write (iunit,'(a)') 'TER'
1207 if (itype(i,1).eq.ntyp1) cycle
1208 if ((itype(i,1).eq.10 .and. itype(i+1,1).ne.ntyp1).or.(molnum(i).eq.5)) then
1209 write (iunit,30) ica(i),ica(i+1)
1210 else if (itype(i,1).ne.10 .and. itype(i+1,1).ne.ntyp1) then
1211 write (iunit,30) ica(i),ica(i+1),ica(i)+1
1212 else if (itype(i,1).ne.10 .and. itype(i+1,1).eq.ntyp1) then
1213 write (iunit,30) ica(i),ica(i)+1
1216 if ((itype(nct,1).ne.10).and.(molnum(i).ne.5)) then
1217 write (iunit,30) ica(nct),ica(nct)+1
1221 write (iunit,30) ica(idssb(i))+1,ica(jdssb(i))+1
1223 write (iunit,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
1226 write (iunit,'(a6)') 'ENDMDL'
1227 10 FORMAT ('ATOM',I7,' CA ',A3,1X,A1,I4,4X,3F8.3,f15.3)
1229 20 FORMAT ('ATOM',I7,' CB ',A3,1X,A1,I4,4X,3F8.3,f15.3)
1230 40 FORMAT ("ATOM",I7," C5' ",1X,2A1,1X,A1,I4,4X,3F8.3,f15.3)
1231 50 FORMAT ("ATOM",I7," C1' ",1X,2A1,1X,A1,I4,4X,3F8.3,f15.3)
1233 30 FORMAT ('CONECT',8I5)
1234 60 FORMAT ('HETATM',I5,' CA ',A3,1X,A1,I4,4X,3F8.3,f15.3)
1237 end subroutine pdbout
1239 !-----------------------------------------------------------------------------
1240 subroutine MOL2out(etot,tytul)
1241 ! Prints the Cartesian coordinates of the alpha-carbons in the Tripos mol2
1243 use geometry_data, only: c
1246 ! implicit real*8 (a-h,o-z)
1247 ! include 'DIMENSIONS'
1248 ! include 'COMMON.CHAIN'
1249 ! include 'COMMON.INTERACT'
1250 ! include 'COMMON.NAMES'
1251 ! include 'COMMON.IOUNITS'
1252 ! include 'COMMON.HEADER'
1253 ! include 'COMMON.SBRIDGE'
1254 character(len=32) :: tytul,fd
1255 character(len=3) :: zahl
1256 character(len=6) :: res_num,pom !,ucase
1260 real(kind=8) :: etot
1264 #elif (defined CRAY)
1269 write (imol2,'(a)') '#'
1270 write (imol2,'(a)') &
1271 '# Creating user name: unres'
1272 write (imol2,'(2a)') '# Creation time: ',&
1274 write (imol2,'(/a)') '\@<TRIPOS>MOLECULE'
1275 write (imol2,'(a)') tytul
1276 write (imol2,'(5i5)') nct-nnt+1,nct-nnt+nss+1,nct-nnt+nss+1,0,0
1277 write (imol2,'(a)') 'SMALL'
1278 write (imol2,'(a)') 'USER_CHARGES'
1279 write (imol2,'(a)') '\@<TRIPOS>ATOM'
1281 write (zahl,'(i3)') i
1282 pom=ucase(restyp(itype(i,1),1))
1283 res_num = pom(:3)//zahl(2:)
1284 write (imol2,10) i,(c(j,i),j=1,3),i,res_num,0.0
1286 write (imol2,'(a)') '\@<TRIPOS>BOND'
1288 write (imol2,'(i5,2i6,i2)') i-nnt+1,i-nnt+1,i-nnt+2,1
1291 write (imol2,'(i5,2i6,i2)') nct-nnt+i,ihpb(i),jhpb(i),1
1293 write (imol2,'(a)') '\@<TRIPOS>SUBSTRUCTURE'
1295 write (zahl,'(i3)') i
1296 pom = ucase(restyp(itype(i,1),1))
1297 res_num = pom(:3)//zahl(2:)
1298 write (imol2,30) i-nnt+1,res_num,i-nnt+1,0
1300 10 FORMAT (I7,' CA ',3F10.4,' C.3',I8,1X,A,F11.4,' ****')
1301 30 FORMAT (I7,1x,A,I14,' RESIDUE',I13,' **** ****')
1303 end subroutine MOL2out
1304 !-----------------------------------------------------------------------------
1308 use energy_data, only: itype
1310 ! implicit real*8 (a-h,o-z)
1311 ! include 'DIMENSIONS'
1312 ! include 'COMMON.IOUNITS'
1313 ! include 'COMMON.CHAIN'
1314 ! include 'COMMON.VAR'
1315 ! include 'COMMON.LOCAL'
1316 ! include 'COMMON.INTERACT'
1317 ! include 'COMMON.NAMES'
1318 ! include 'COMMON.GEO'
1319 ! include 'COMMON.TORSION'
1323 write (iout,'(/a)') 'Geometry of the virtual chain.'
1324 write (iout,'(7a)') ' Res ',' d',' Theta',&
1325 ' Phi',' Dsc',' Alpha',' Omega'
1328 ! print *,vbld(i),"vbld(i)",i
1329 write (iout,'(a3,i6,6f10.3)') restyp(iti,1),i,vbld(i),&
1330 rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),rad2deg*alph(i),&
1334 end subroutine intout
1335 !-----------------------------------------------------------------------------
1337 subroutine briefout(it,ener,free)!,plik)
1339 subroutine briefout(it,ener)
1344 ! implicit real*8 (a-h,o-z)
1345 ! include 'DIMENSIONS'
1346 ! include 'COMMON.IOUNITS'
1347 ! include 'COMMON.CHAIN'
1348 ! include 'COMMON.VAR'
1349 ! include 'COMMON.LOCAL'
1350 ! include 'COMMON.INTERACT'
1351 ! include 'COMMON.NAMES'
1352 ! include 'COMMON.GEO'
1353 ! include 'COMMON.SBRIDGE'
1354 ! print '(a,i5)',intname,igeom
1357 real(kind=8) :: ener,free
1358 ! character(len=80) :: plik
1361 #if defined(AIX) || defined(PGI)
1362 open (igeom,file=intname,position='append')
1364 open (igeom,file=intname,access='append')
1373 WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,NSS)
1375 WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,8)
1377 WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,NSS)
1379 WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,9)
1381 WRITE (igeom,190) (IHPB(I),JHPB(I),I=10,NSS)
1383 ! IF (nvar.gt.nphi) WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
1384 WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
1385 WRITE (igeom,200) (RAD2DEG*PHI(I),I=4,NRES)
1386 ! if (nvar.gt.nphi+ntheta) then
1387 write (igeom,200) (rad2deg*alph(i),i=2,nres-1)
1388 write (igeom,200) (rad2deg*omeg(i),i=2,nres-1)
1391 180 format (I5,F12.3,I2,9(1X,2I3))
1392 190 format (3X,11(1X,2I3))
1395 end subroutine briefout
1396 !-----------------------------------------------------------------------------
1398 subroutine fdate(fd)
1399 character(len=32) :: fd
1402 end subroutine fdate
1404 !-----------------------------------------------------------------------------
1406 real(kind=8) function gyrate(jcon)
1408 real(kind=8) function gyrate()
1411 use geometry_data, only: c
1414 ! implicit real*8 (a-h,o-z)
1415 ! include 'DIMENSIONS'
1416 ! include 'COMMON.INTERACT'
1417 ! include 'COMMON.CHAIN'
1419 real(kind=8),dimension(3) :: cen
1429 cen(j)=cen(j)+c(j,i)
1433 cen(j)=cen(j)/dble(nct-nnt+1)
1438 rg = rg + (c(j,i)-cen(j))**2
1442 gyrate = dsqrt(rg/dble(nct-nnt+1))
1444 gyrate = sqrt(rg/dble(nct-nnt+1))
1449 !-----------------------------------------------------------------------------
1451 subroutine reads(rekord,lancuch,wartosc,default)
1453 character*(*) :: rekord,lancuch,wartosc,default
1454 character(len=80) :: aux
1455 integer :: lenlan,lenrec,iread,ireade
1459 lenlan=ilen(lancuch)
1461 iread=index(rekord,lancuch(:lenlan)//"=")
1462 ! print *,"rekord",rekord," lancuch",lancuch
1463 ! print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
1464 if (iread.eq.0) then
1468 iread=iread+lenlan+1
1469 ! print *,"iread",iread
1470 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
1471 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
1473 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
1475 ! print *,"iread",iread
1476 if (iread.gt.lenrec) then
1481 ! print *,"ireade",ireade
1482 do while (ireade.lt.lenrec .and. &
1483 .not.iblnk(rekord(ireade:ireade)))
1486 wartosc=rekord(iread:ireade)
1488 end subroutine reads
1490 !-----------------------------------------------------------------------------
1492 !-----------------------------------------------------------------------------
1493 subroutine permut(isym,nperm,tabperm)
1494 !c integer maxperm,maxsym
1495 !c parameter (maxperm=3628800)
1496 !c parameter (maxsym=10)
1497 ! include "DIMENSIONS"
1498 integer n,a,tabperm,nperm,kkk,i,isym
1501 dimension a(isym),tabperm(50,5040)
1516 !c print '(i3,2x,100i3)',kkk+1,(a(i),i=1,n)
1521 if(nextp(n,a)) go to 10
1526 !-----------------------------------------------------------------------------
1527 logical function nextp(n,a)
1529 integer :: n,i,j,k,t
1531 integer,dimension(n) :: a
1533 10 if(a(i).lt.a(i+1)) go to 20
1550 if(a(j).lt.a(i)) go to 40
1558 !-----------------------------------------------------------------------------
1560 !-----------------------------------------------------------------------------
1561 integer function rescode(iseq,nam,itype,molecule)
1563 ! use io_base, only: ucase
1564 ! implicit real*8 (a-h,o-z)
1565 ! include 'DIMENSIONS'
1566 ! include 'COMMON.NAMES'
1567 ! include 'COMMON.IOUNITS'
1568 character(len=3) :: nam !,ucase
1569 integer :: iseq,itype,i
1571 print *,molecule,nam
1572 if (molecule.eq.1) then
1573 if (itype.eq.0) then
1575 do i=-ntyp1_molec(molecule),ntyp1_molec(molecule)
1576 if (ucase(nam).eq.restyp(i,molecule)) then
1584 do i=-ntyp1_molec(molecule),ntyp1_molec(molecule)
1585 if (nam(1:1).eq.onelet(i)) then
1592 else if (molecule.eq.2) then
1593 do i=1,ntyp1_molec(molecule)
1594 print *,nam(1:1),restyp(i,molecule)(1:1)
1595 print *,nam(2:2),"read",i
1596 if (nam(2:2).eq.restyp(i,molecule)(1:1)) then
1601 else if (molecule.eq.3) then
1602 write(iout,*) "SUGAR not yet implemented"
1604 else if (molecule.eq.4) then
1605 write(iout,*) "Explicit LIPID not yet implemented"
1607 else if (molecule.eq.5) then
1608 do i=-1,ntyp1_molec(molecule)
1609 print *,restyp(i,molecule)
1610 print *,i,restyp(i,molecule)(1:2),nam(1:2)
1611 if (ucase(nam(1:2)).eq.restyp(i,molecule)(1:2)) then
1617 write(iout,*) "molecule not defined"
1619 write (iout,10) iseq,nam
1621 10 format ('**** Error - residue',i4,' has an unresolved name ',a3)
1622 end function rescode
1623 integer function sugarcode(sugar,ires)
1626 if (sugar.eq.'D') then
1628 else if (sugar.eq.' ') then
1631 write (iout,*) 'UNKNOWN sugar type for residue',ires,' ',sugar
1635 end function sugarcode
1636 integer function rescode_lip(res,ires)
1637 character(len=3):: res
1640 if (res.eq.'NC3') rescode_lip =4
1641 if (res.eq.'NH3') rescode_lip =2
1642 if (res.eq.'PO4') rescode_lip =3
1643 if (res.eq.'GL1') rescode_lip =12
1644 if (res.eq.'GL2') rescode_lip =12
1645 if (res.eq.'C1A') rescode_lip =18
1646 if (res.eq.'C2A') rescode_lip =18
1647 if (res.eq.'C3A') rescode_lip =18
1648 if (res.eq.'C4A') rescode_lip =18
1649 if (res.eq.'C1B') rescode_lip =18
1650 if (res.eq.'C2B') rescode_lip =18
1651 if (res.eq.'C3B') rescode_lip =18
1652 if (res.eq.'C4B') rescode_lip =18
1653 if (res.eq.'D2A') rescode_lip =16
1654 if (res.eq.'D2B') rescode_lip =16
1656 if (rescode_lip.eq.0) write(iout,*) "UNKNOWN RESIDUE",ires,res
1659 !-----------------------------------------------------------------------------
1660 !-----------------------------------------------------------------------------