2 !-----------------------------------------------------------------------------
5 use geometry, only:int_from_cart,sc_loc_geom
6 use io_config,only:readpdb,readpdb_template
8 use io_base !, only: ilen
9 use geometry_data, only: nres,c,nres_molec,boxxsize,boxysize,boxzsize
10 use energy_data, only: nnt,nct,nss
11 use control_data, only: lside
13 !-----------------------------------------------------------------------------
16 !-----------------------------------------------------------------------------
18 !-----------------------------------------------------------------------------
20 !-----------------------------------------------------------------------------
21 SUBROUTINE WRTCLUST(NCON,ICUT,PRINTANG,PRINTPDB,printmol2,ib)
23 use hc_, only: ioffset
24 use control_data, only: lprint_cart,lprint_int,titel
25 use geometry, only: int_from_cart1,nres
26 ! implicit real*8 (a-h,o-z)
27 ! include 'DIMENSIONS'
28 ! include 'sizesclu.dat'
29 integer,parameter :: num_in_line=5
30 LOGICAL :: PRINTANG(max_cut)
31 integer :: PRINTPDB(max_cut),printmol2(max_cut)
32 ! include 'COMMON.CONTROL'
33 ! include 'COMMON.HEADER'
34 ! include 'COMMON.CHAIN'
35 ! include 'COMMON.VAR'
36 ! include 'COMMON.CLUSTER'
37 ! include 'COMMON.IOUNITS'
38 ! include 'COMMON.GEO'
39 ! include 'COMMON.FREE'
40 ! include 'COMMON.TEMPFAC'
41 real(kind=8) :: rmsave(maxgr)
42 CHARACTER(len=64) :: prefixp,NUMM,MUMM,EXTEN,extmol
43 character(len=80) :: cfname
44 character(len=8) :: ctemper
45 DATA EXTEN /'.pdb'/,extmol /'.mol2'/,NUMM /'000'/,&
48 integer :: ncon,icut,ib
49 integer :: i,ii,ii1,ii2,igr,ind1,ind2,ico,icon,&
50 irecord,nrecord,j,k,jj,ind,ncon_lim,ncon_out
51 real(kind=8) :: temper,curr_dist,emin,qpart,boltz,&
52 ave_dim,amax_dim,emin1
54 if (allocated(tempfac)) deallocate(tempfac)
55 allocate(tempfac(2,nres))
60 ! print *,"calling WRTCLUST",ncon
61 ! write (iout,*) "ICUT",icut," PRINTPDB ",PRINTPDB(icut)
64 temper=1.0d0/(beta_h(ib)*1.987d-3)
65 if (temper.lt.100.0d0) then
66 write(ctemper,'(f3.0)') temper
68 else if (temper.lt.1000.0) then
69 write (ctemper,'(f4.0)') temper
72 write (ctemper,'(f5.0)') temper
76 do i=1,ncon*(ncon-1)/2
79 close(80,status='delete')
81 ! PRINT OUT THE RESULTS OF CLUSTER ANALYSIS
83 ii1= index(intinname,'/')
88 ii2=index(intinname(ii1:),'/')
90 ii = ii1+index(intinname(ii1:),'.')-1
96 prefixp=intinname(ii1:ii)
97 !d print *,icut,printang(icut),printpdb(icut),printmol2(icut)
98 !d print *,'ecut=',ecut
101 WRITE (iout,200) IGR,totfree_gr(igr)/beta_h(ib),LICZ(IGR)
102 NRECORD=LICZ(IGR)/num_in_line
104 DO 63 IRECORD=1,NRECORD
105 IND2=IND1+num_in_line-1
106 WRITE (iout,300) (list_conf(NCONF(IGR,ICO)),&
107 totfree(NCONF(IGR,ICO))/beta_h(ib),ICO=IND1,IND2)
110 WRITE (iout,300) (list_conf(NCONF(IGR,ICO)),&
111 totfree(NCONF(IGR,ICO))/beta_h(ib),ICO=IND1,LICZ(IGR))
113 ICON=list_conf(NCONF(IGR,1))
114 ! WRITE (iout,'(16F5.0)') (rad2deg*phiall(IND,icon),IND=4,nphi+3)
115 ! 12/8/93 Estimation of "diameters" of the subsequent families.
118 ! write (iout,*) "ecut",ecut
121 if (totfree(ii)-emin .gt. ecut) goto 10
126 ind=ioffset(ncon,ii,jj)
128 ind=ioffset(ncon,jj,ii)
130 ! write (iout,*) " ncon",ncon,"i",i," j",j," ii",ii," jj",jj,
133 curr_dist=dabs(diss(ind)+0.0d0)
134 ! write(iout,'(i10,4i4,f12.4)') ind,ii,jj,list_conf(ii),
135 ! & list_conf(jj),curr_dist
136 if (curr_dist .gt. amax_dim) amax_dim=curr_dist
137 ave_dim=ave_dim+curr_dist**2
140 10 if (licz(igr) .gt. 1) &
141 ave_dim=sqrt(ave_dim/(licz(igr)*(licz(igr)-1)/2))
142 write (iout,'(/A,F8.1,A,F8.1)') &
143 'Max. distance in the family:',amax_dim,&
144 '; average distance in the family:',ave_dim
149 boltz=dexp(-totfree(icon))
150 rmsave(igr)=rmsave(igr)+boltz*rmstb(icon)
153 rmsave(igr)=rmsave(igr)/qpart
154 write (iout,'(a,f5.2,a)') "Average RMSD",rmsave(igr)," A"
157 WRITE (iout,500) (list_conf(I),IASS(I),I=1,NCON)
158 ! print *,icut,printang(icut)
159 IF (PRINTANG(ICUT) .and. (lprint_cart .or. lprint_int)) then
161 ! print *,'emin',emin,' ngr',ngr
162 if (lprint_cart) then
163 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))&
166 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))&
171 if (totfree_gr(igr)-emin.le.ecut) then
172 if (lprint_cart) then
173 call cartout(igr,icon,totfree(icon)/beta_h(ib),&
174 totfree_gr(igr)/beta_h(ib),&
177 ! print '(a)','calling briefout'
180 c(j,i)=allcart(j,i,icon)
183 call int_from_cart1(.false.)
184 !el call briefout(igr,iscore(icon),totfree(icon)/beta_h(ib),&
185 !el totfree_gr(igr),nss_all(icon),ihpb_all(1,icon),&
186 !el jhpb_all(1,icon),cfname)
187 call briefout(igr,totfree(icon)/beta_h(ib),&
189 ! print '(a)','exit briefout'
195 IF (PRINTPDB(ICUT).gt.0) THEN
196 ! Write out a number of conformations from each family in PDB format and
197 ! create InsightII command file for their displaying in different colors
198 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))&
200 write (iout,*) "cfname",cfname
201 OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
202 write (ipdb,'(a,f8.2)') &
203 "REMAR AVERAGE CONFORMATIONS AT TEMPERATURE",temper
209 DO WHILE(I.LE.NGR .AND. totfree_gr(i)-EMIN.LE.ECUT)
210 ! write (iout,*) "i",i," ngr",ngr,totfree_gr(I),EMIN,ecut
211 write (NUMM,'(bz,i4.4)') i
212 ncon_lim=min0(licz(i),printpdb(icut))
213 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))&
214 //"K_"//numm(:ilen(numm))//exten
215 OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
216 write (ipdb,'("REMARK CLUSTER",i5," FREE ENERGY",1pe14.5," AVE RMSD",0pf5.2)')&
217 i,totfree_gr(i)/beta_h(ib),rmsave(i) !'
218 ! Write conformations of the family i to PDB files
220 do while (ncon_out.lt.printpdb(icut) .and. &
221 ncon_out.lt.licz(i).and. &
222 totfree(nconf(i,ncon_out+1))-EMIN1.LE.ECUT)
224 ! write (iout,*) i,ncon_out,nconf(i,ncon_out),
225 ! & totfree(nconf(i,ncon_out)),emin1,ecut
227 write (iout,*) "ncon_out",ncon_out
237 c(k,ii)=allcart(k,ii,icon)
240 call pdboutC(totfree(icon)/beta_h(ib),rmstb(icon),titel)
241 write (ipdb,'("TER")')
244 ! Average structures and structures closest to average
245 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))&
247 OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED',&
250 write (ipdb,'(a,i5)') "REMARK CLUSTER",i
251 call pdboutC(totfree_gr(i)/beta_h(ib),rmsave(i),titel)
252 write (ipdb,'("TER")')
253 call closest_coord(i)
254 call pdboutC(totfree_gr(i)/beta_h(ib),rmsave(i),titel)
255 write (ipdb,'("TER")')
262 IF (printmol2(icut).gt.0) THEN
263 ! Write out a number of conformations from each family in PDB format and
264 ! create InsightII command file for their displaying in different colors
269 DO WHILE(I.LE.NGR .AND. totfree_gr(i)-EMIN.LE.ECUT)
270 write (NUMM,'(bz,i4.4)') i
271 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))&
272 //"K_"//numm(:ilen(numm))//extmol
273 OPEN(imol2,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
275 do while (ncon_out.lt.printmol2(icut) .and. &
276 ncon_out.lt.licz(i).and. &
277 totfree(nconf(i,ncon_out+1))-EMIN1.LE.ECUT)
284 c(k,ii)=allcart(k,ii,icon)
287 CALL MOL2OUT(totfree(icon)/beta_h(ib),'STRUCTURE'//numm)
295 100 FORMAT (//'THERE ARE ',I4,' FAMILIES OF CONFORMATIONS')
296 200 FORMAT (/'FAMILY ',I4,' WITH TOTAL FREE ENERGY',1pE15.5,&
297 ' CONTAINS ',I4,' CONFORMATION(S): ')
298 ! 300 FORMAT ( 8(I4,F6.1))
299 300 FORMAT (5(I4,1pe12.3))
300 400 FORMAT (//'ASSIGNMENT OF CONSECUTIVE CONFORMATIONS TO FAMILIES:')
301 500 FORMAT (8(2I4,2X))
302 600 FORMAT ('REMARK FAMILY',I4,' CONFORMATION',I4,' ENERGY ',E15.6)
304 END SUBROUTINE WRTCLUST
305 !------------------------------------------------------------------------------
306 !------------------------------------------------------------------------------
307 subroutine ave_coord(igr)
309 use control_data, only:lside
310 use regularize_, only:fitsq,matvec
312 ! include 'DIMENSIONS'
313 ! include 'sizesclu.dat'
314 ! include 'COMMON.CONTROL'
315 ! include 'COMMON.CLUSTER'
316 ! include 'COMMON.CHAIN'
317 ! include 'COMMON.INTERACT'
318 ! include 'COMMON.VAR'
319 ! include 'COMMON.TEMPFAC'
320 ! include 'COMMON.IOUNITS'
322 real(kind=8) :: przes(3),obrot(3,3)
323 real(kind=8) :: xx(3,2*nres),yy(3,2*nres),csq(3,2*nres) !(3,maxres2)
325 integer :: i,ii,j,k,icon,jcon,igr
326 real(kind=8) :: rms,boltz,qpart,cwork(3,2*nres),cref1(3,2*nres) !(3,maxres2)
327 ! write (iout,*) "AVE_COORD: igr",igr
330 boltz = dexp(-totfree(jcon)+eref)
334 c(j,i)=allcart(j,i,jcon)*boltz
335 cref1(j,i)=allcart(j,i,jcon)
336 csq(j,i)=allcart(j,i,jcon)**2*boltz
346 xx(j,ii)=allcart(j,i,jcon)
351 ! if (itype(i).ne.10) then
354 xx(j,ii)=allcart(j,i+nres,jcon)
355 yy(j,ii)=cref1(j,i+nres)
359 call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
363 cwork(j,i)=allcart(j,i,jcon)
366 call fitsq(rms,cwork(1,nnt),cref1(1,nnt),nct-nnt+1,przes,obrot &
369 ! write (iout,*) "rms",rms
371 ! write (iout,'(i3,f10.5,5x,3f10.5)')i,przes(i),(obrot(i,j),j=1,3)
374 print *,'error, rms^2 = ',rms,icon,jcon
377 if (non_conv) print *,non_conv,icon,jcon
378 boltz=dexp(-totfree(jcon)+eref)
379 qpart = qpart + boltz
382 xx(j,i)=allcart(j,i,jcon)
385 call matvec(cwork,obrot,xx,2*nres)
387 ! write (iout,'(i5,2(3f10.5,5x))') i,(cwork(j,i),j=1,3),
388 ! & (allcart(j,i,jcon),j=1,3)
390 cwork(j,i)=cwork(j,i)+przes(j)
391 c(j,i)=c(j,i)+cwork(j,i)*boltz
392 csq(j,i)=csq(j,i)+cwork(j,i)**2*boltz
399 csq(j,i)=csq(j,i)/qpart-c(j,i)**2
401 ! write (iout,'(i5,3f10.5)') i,(csq(j,i),j=1,3)
407 tempfac(1,i)=tempfac(1,i)+csq(j,i)
408 tempfac(2,i)=tempfac(2,i)+csq(j,i+nres)
410 tempfac(1,i)=dsqrt(tempfac(1,i))
411 tempfac(2,i)=dsqrt(tempfac(2,i))
414 end subroutine ave_coord
415 !------------------------------------------------------------------------------
416 subroutine closest_coord(igr)
418 use regularize_, only:fitsq
420 ! include 'DIMENSIONS'
421 ! include 'sizesclu.dat'
422 ! include 'COMMON.IOUNITS'
423 ! include 'COMMON.CONTROL'
424 ! include 'COMMON.CLUSTER'
425 ! include 'COMMON.CHAIN'
426 ! include 'COMMON.INTERACT'
427 ! include 'COMMON.VAR'
429 real(kind=8) :: przes(3),obrot(3,3)
430 real(kind=8) :: xx(3,2*nres),yy(3,2*nres) !(3,maxres2)
431 integer :: i,ii,j,k,icon,jcon,jconmin,igr
432 real(kind=8) :: rms,rmsmin,cwork(3,2*nres)
442 xx(j,ii)=allcart(j,i,jcon)
447 ! if (itype(i).ne.10) then
450 xx(j,ii)=allcart(j,i+nres,jcon)
455 call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
459 cwork(j,i)=allcart(j,i,jcon)
462 call fitsq(rms,cwork(1,nnt),c(1,nnt),nct-nnt+1,przes,obrot,&
466 print *,'error, rms^2 = ',rms,icon,jcon
469 ! write (iout,*) "jcon",jcon," rms",rms," rmsmin",rmsmin
470 if (non_conv) print *,non_conv,icon,jcon
471 if (rms.lt.rmsmin) then
476 ! write (iout,*) "rmsmin",rmsmin," rms",rms
480 c(j,i)=allcart(j,i,jconmin)
484 end subroutine closest_coord
485 !-----------------------------------------------------------------------------
487 !-----------------------------------------------------------------------------
488 subroutine read_coords(ncon,*)
490 use energy_data, only: ihpb,jhpb,max_ene
491 use control_data, only: from_bx,from_cx
492 use control, only: tcpu
494 ! include "DIMENSIONS"
495 ! include "sizesclu.dat"
499 integer :: IERROR,ERRCODE !,STATUS(MPI_STATUS_SIZE)
500 ! include "COMMON.MPI"
502 use MPI_data, only: nprocs
504 ! include "COMMON.CONTROL"
505 ! include "COMMON.CHAIN"
506 ! include "COMMON.INTERACT"
507 ! include "COMMON.IOUNITS"
508 ! include "COMMON.VAR"
509 ! include "COMMON.SBRIDGE"
510 ! include "COMMON.GEO"
511 ! include "COMMON.CLUSTER"
512 character(len=3) :: liczba
514 integer :: i,j,jj,jjj,jj_old,icount,k,kk,l,ii,if,ib,&
516 integer :: ixdrf,iret,itmp
517 real(kind=4) :: prec,reini,refree,rmsdev
518 integer :: nrec,nlines,iscor,lenrec,lenrec_in
519 real(kind=8) :: energ,t_acq !,tcpu
520 !el integer ilen,iroof
521 !el external ilen,iroof
522 real(kind=8) :: rjunk
523 integer :: ntot_all(0:nprocs-1) !(0:maxprocs-1)
525 real(kind=8) :: energia(0:max_ene),etot
526 real(kind=4) :: csingle(3,2*nres+2)
527 integer :: Previous,Next
528 character(len=256) :: bprotfiles
529 ! print *,"Processor",me," calls read_protein_data"
531 if (me.eq.master) then
532 Previous=MPI_PROC_NULL
536 if (me.eq.nprocs-1) then
541 ! Set the scratchfile names
542 write (liczba,'(bz,i3.3)') me
544 allocate(STATUS(MPI_STATUS_SIZE))
546 ! 1/27/05 AL Change stored coordinates to single precision and don't store
547 ! energy components in the binary databases.
548 lenrec=12*(nres+nct-nnt+1)+4*(2*nss+2)+16
549 lenrec_in=12*(nres+nct-nnt+1)+4*(2*nss+2)+24
551 write (iout,*) "nres",nres," nnt",nnt," nct",nct," nss", nss
552 write (iout,*) "lenrec_in",lenrec_in
554 bprotfiles=scratchdir(:ilen(scratchdir))// &
555 "/"//prefix(:ilen(prefix))//liczba//".xbin"
557 ! allocate cluster arrays
558 allocate(energy(0:maxconf),totfree(0:maxconf)) !(0:maxconf)
559 allocate(entfac(maxconf)) !(maxconf)
560 allocate(rmstb(maxconf)) !(maxconf)
561 allocate(allcart(3,2*nres,maxstr_proc)) !(3,maxres2,maxstr_proc)
562 allocate(nss_all(maxstr_proc)) !(maxstr_proc)
563 allocate(ihpb_all(maxss,maxstr_proc),jhpb_all(maxss,maxstr_proc))!(maxss,maxstr_proc)
564 allocate(iscore(maxconf)) !(maxconf)
572 if (from_cart .and. .not. from_bx .and. .not. from_cx) then
574 read (intin,*,end=13,err=11) energy(icon),totfree(icon),&
576 nss_all(icon),(ihpb_all(ii,icon),jhpb_all(i,icon),&
577 i=1,nss_all(icon)),iscore(icon)
579 read (intin,*,end=13,err=11) energy(icon),rmstb(icon),&
580 nss_all(icon),(ihpb_all(ii,icon),jhpb_all(i,icon),&
581 i=1,nss_all(icon)),iscore(icon)
583 read (intin,'(8f10.5)',end=13,err=10) &
584 ((allcart(j,i,icon),j=1,3),i=1,nres),&
585 ((allcart(j,i+nres,icon),j=1,3),i=nnt,nct)
586 print *,icon,energy(icon),nss_all(icon),rmstb(icon)
588 read(intin,'(a80)',end=13,err=12) lineh
589 read(lineh(:5),*,err=8) ic
591 read(lineh(6:),*,err=8) energy(icon)
593 read(lineh(6:),*,err=8) energy(icon)
597 print *,'error, assuming e=1d10',lineh
601 !old read(lineh(18:),*,end=13,err=11) nss_all(icon)
602 ii = index(lineh(15:)," ")+15
603 read(lineh(ii:),*,end=13,err=11) nss_all(icon)
604 IF (NSS_all(icon).LT.9) THEN
605 read (lineh(20:),*,end=102) &
606 (IHPB_all(I,icon),JHPB_all(I,icon),I=1,NSS_all(icon)),&
609 read (lineh(20:),*,end=102) &
610 (IHPB_all(I,icon),JHPB_all(I,icon),I=1,8)
611 read (intin,*) (IHPB_all(I,icon),JHPB_all(I,icon),&
612 I=9,NSS_all(icon)),iscore(icon)
617 PRINT *,'IC:',IC,' ENERGY:',ENERGY(ICON)
618 call read_angles(intin,*13)
620 phiall(i,icon)=phi(i)
621 thetall(i,icon)=theta(i)
622 alphall(i,icon)=alph(i)
623 omall(i,icon)=omeg(i)
629 ! CALCULATE DISTANCES
631 10 print *,'something wrong with angles'
633 11 print *,'something wrong with NSS',nss
635 12 print *,'something wrong with header'
642 open (icbase,file=bprotfiles,status="unknown",&
643 form="unformatted",access="direct",recl=lenrec)
644 ! Read conformations from binary DA files (one per batch) and write them to
645 ! a binary DA scratchfile.
649 write (liczba,'(bz,i3.3)') me
650 IF (ME.EQ.MASTER) THEN
651 ! Only the master reads the database; it'll send it to the other procs
659 open (intin,file=intinname,status="old",form="unformatted",&
660 access="direct",recl=lenrec_in)
662 else if (from_cx) then
663 #if (defined(AIX) && !defined(JUBL))
664 call xdrfopen_(ixdrf,intinname, "r", iret)
666 call xdrfopen(ixdrf,intinname, "r", iret)
669 write (iout,*) "xdrfopen: iret",iret
671 write (iout,*) "Error: coordinate file ",&
672 intinname(:ilen(intinname))," does not exist."
675 call MPI_ABORT(MPI_COMM_WORLD,IERROR,ERRCODE)
680 write (iout,*) "Error: coordinate format not specified"
683 call MPI_ABORT(MPI_COMM_WORLD,IERROR,ERRCODE)
691 write (iout,*) "Opening file ",intinname(:ilen(intinname))
692 write (iout,*) "lenrec",lenrec_in
696 ! write (iout,*) "maxconf",maxconf
700 !el if (i.gt.maxconf) then
701 !el write (iout,*) "Error: too many conformations ",&
702 !el "(",maxconf,") maximum."
704 !el call MPI_Abort(MPI_COMM_WORLD,errcode,ierror)
708 ! write (iout,*) "i",i
711 read(intin,err=101,end=101) &
712 ((csingle(l,k),l=1,3),k=1,nres),&
713 ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
714 nss,(ihpb(k),jhpb(k),k=1,nss),&
716 entfac(jj+1),rmstb(jj+1),iscor
723 #if (defined(AIX) && !defined(JUBL))
724 call xdrf3dfcoord_(ixdrf, csingle, itmp, prec, iret)
725 if (iret.eq.0) goto 101
726 call xdrfint_(ixdrf, nss, iret)
727 if (iret.eq.0) goto 101
729 call xdrfint_(ixdrf, ihpb(j), iret)
730 if (iret.eq.0) goto 101
731 call xdrfint_(ixdrf, jhpb(j), iret)
732 if (iret.eq.0) goto 101
734 call xdrffloat_(ixdrf,reini,iret)
735 if (iret.eq.0) goto 101
736 call xdrffloat_(ixdrf,refree,iret)
737 if (iret.eq.0) goto 101
738 call xdrffloat_(ixdrf,rmsdev,iret)
739 if (iret.eq.0) goto 101
740 call xdrfint_(ixdrf,iscor,iret)
741 if (iret.eq.0) goto 101
743 ! write (iout,*) "calling xdrf3dfcoord"
744 call xdrf3dfcoord(ixdrf, csingle, itmp, prec, iret)
745 ! write (iout,*) "iret",iret
747 if (iret.eq.0) goto 101
748 call xdrfint(ixdrf, nss, iret)
749 write (iout,*) "iret",iret
750 ! write (iout,*) "nss",nss,i,"TUTU"
752 if (iret.eq.0) goto 101
754 call xdrfint(ixdrf, ihpb(k), iret)
755 if (iret.eq.0) goto 101
756 call xdrfint(ixdrf, jhpb(k), iret)
757 if (iret.eq.0) goto 101
758 ! write(iout,*), ihpb(k),jhpb(k),"TUTU"
760 call xdrffloat(ixdrf,reini,iret)
761 if (iret.eq.0) goto 101
762 ! write(iout,*) "TUTU", reini
763 call xdrffloat(ixdrf,refree,iret)
764 ! write(iout,*) "TUTU", refree
765 if (iret.eq.0) goto 101
766 call xdrffloat(ixdrf,rmsdev,iret)
767 if (iret.eq.0) goto 101
768 call xdrfint(ixdrf,iscor,iret)
769 if (iret.eq.0) goto 101
781 c(l,nres+k)=csingle(l,nres+k-nnt+1)
786 write (iout,'(5hREAD ,i5,3f15.4,i10)') &
787 jj+1,energy(jj+1),entfac(jj+1),&
789 write (iout,*) "Conformation",jjj+1,jj+1
790 write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
791 write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
794 call add_new_cconf(jjj,jj,jj_old,icount,Next)
797 write (iout,*) i-1," conformations read from DA file ",&
798 intinname(:ilen(intinname))
799 write (iout,*) jj," conformations read so far"
803 #if (defined(AIX) && !defined(JUBL))
804 call xdrfclose_(ixdrf, iret)
806 call xdrfclose(ixdrf, iret)
811 write (iout,*) "jj_old",jj_old," jj",jj
813 call write_and_send_cconf(icount,jj_old,jj,Next)
814 call MPI_Send(0,1,MPI_INTEGER,Next,570,&
815 MPI_COMM_WORLD,IERROR)
818 call write_and_send_cconf(icount,jj_old,jj,Next)
820 t_acq = tcpu() - t_acq
822 write (iout,*) "Processor",me,&
823 " time for conformation read/send",t_acq
825 ! A worker gets the confs from the master and sends them to its neighbor
827 call receive_and_pass_cconf(icount,jj_old,jj,&
829 t_acq = tcpu() - t_acq
836 write(iout,*)"A total of",ncon," conformations read."
838 allocate(enetb(1:max_ene,ncon)) !(max_ene,maxstr_proc)
840 ! Check if everyone has the same number of conformations
841 call MPI_Allgather(ncon,1,MPI_INTEGER,&
842 ntot_all(0),1,MPI_INTEGER,MPI_Comm_World,IERROR)
846 if (ncon.ne.ntot_all(i)) then
847 write (iout,*) "Number of conformations at processor",i,&
848 " differs from that at processor",me,&
856 write (iout,*) "Number of conformations read by processors"
859 write (iout,'(8i10)') i,ntot_all(i)
861 write (iout,*) "Calculation terminated."
867 1111 write(iout,*) "Error opening coordinate file ",&
868 intinname(:ilen(intinname))
871 end subroutine read_coords
872 !------------------------------------------------------------------------------
873 subroutine add_new_cconf(jjj,jj,jj_old,icount,Next)
875 use geometry_data, only: vbld,rad2deg,theta,phi,alph,omeg,deg2rad
876 use energy_data, only: itel,itype,dsc,max_ene,molnum
877 use control_data, only: symetr
878 use geometry, only: int_from_cart1
880 ! include "DIMENSIONS"
881 ! include "sizesclu.dat"
882 ! include "COMMON.CLUSTER"
883 ! include "COMMON.CONTROL"
884 ! include "COMMON.CHAIN"
885 ! include "COMMON.INTERACT"
886 ! include "COMMON.LOCAL"
887 ! include "COMMON.IOUNITS"
888 ! include "COMMON.NAMES"
889 ! include "COMMON.VAR"
890 ! include "COMMON.SBRIDGE"
891 ! include "COMMON.GEO"
892 integer :: i,j,jj,jjj,jj_old,icount,k,kk,l,ii,ib,&
893 nn,nn1,inan,Next,itj,chalen,mnum
894 real(kind=8) :: etot,energia(0:max_ene)
896 chalen=int((nct-nnt+2)/symetr)
897 call int_from_cart1(.false.)
900 write (iout,*) "Check atom",j,itype(j,mnum),vbld(j)
902 if (itype(j,mnum).eq.ntyp1_molec(mnum)) cycle
903 if (itype(j-1,molnum(j-1)).eq.ntyp1_molec(molnum(j-1))) cycle
905 if ((vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0).and.(mnum.ne.5)) then
907 if (itel(j).ne.0 .and. itel(j-1).ne.0) then
908 write (iout,*) "Conformation",jjj,jj+1,itype(j-1,molnum(j-1)),itype(j,mnum)
909 write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),itel(j),&
911 write (iout,*) "The Cartesian geometry is:"
912 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
913 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
914 write (iout,*) "The internal geometry is:"
915 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
916 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
917 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
918 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
919 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
920 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
922 "This conformation WILL NOT be added to the database."
932 if (itype(j,1).ne.10 .and. (itype(j,mnum).ne.ntyp1_molec(mnum)) &
933 .and. (vbld(nres+j)-dsc(iabs(itj))) &
935 write (iout,*) "Conformation",jjj,jj+1
936 write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j)
937 write (iout,*) "The Cartesian geometry is:"
938 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
939 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
940 write (iout,*) "The internal geometry is:"
941 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
942 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
943 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
944 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
945 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
946 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
948 "This conformation WILL NOT be added to the database."
953 if (theta(j).le.0.0d0) then
955 "Zero theta angle(s) in conformation",jjj,jj+1
956 write (iout,*) "The Cartesian geometry is:"
957 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
958 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
959 write (iout,*) "The internal geometry is:"
960 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
961 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
962 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
963 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
964 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
965 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
967 "This conformation WILL NOT be added to the database."
970 if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad
974 write (iout,*) "Conformation",jjj,jj
975 write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
976 write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
977 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
978 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
979 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
980 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
981 write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct)
982 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
983 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
984 write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
985 write (iout,'(e15.5,16i5)') entfac(icount+1)
986 ! & iscore(icount+1,0)
989 call store_cconf_from_file(jj,icount)
990 if (icount.eq.maxstr_proc) then
992 write (iout,* ) "jj_old",jj_old," jj",jj
994 call write_and_send_cconf(icount,jj_old,jj,Next)
999 end subroutine add_new_cconf
1000 !------------------------------------------------------------------------------
1001 subroutine store_cconf_from_file(jj,icount)
1003 use energy_data, only: ihpb,jhpb
1005 ! include "DIMENSIONS"
1006 ! include "sizesclu.dat"
1007 ! include "COMMON.CLUSTER"
1008 ! include "COMMON.CHAIN"
1009 ! include "COMMON.SBRIDGE"
1010 ! include "COMMON.INTERACT"
1011 ! include "COMMON.IOUNITS"
1012 ! include "COMMON.VAR"
1013 integer :: i,j,jj,icount
1014 ! Store the conformation that has been read in
1017 allcart(j,i,icount)=c(j,i)
1022 ihpb_all(i,icount)=ihpb(i)
1023 jhpb_all(i,icount)=jhpb(i)
1026 end subroutine store_cconf_from_file
1027 !------------------------------------------------------------------------------
1028 subroutine write_and_send_cconf(icount,jj_old,jj,Next)
1031 ! include "DIMENSIONS"
1032 ! include "sizesclu.dat"
1037 ! include "COMMON.MPI"
1039 ! include "COMMON.CHAIN"
1040 ! include "COMMON.SBRIDGE"
1041 ! include "COMMON.INTERACT"
1042 ! include "COMMON.IOUNITS"
1043 ! include "COMMON.CLUSTER"
1044 ! include "COMMON.VAR"
1045 integer :: icount,jj_old,jj,Next
1046 ! Write the structures to a scratch file
1048 ! Master sends the portion of conformations that have been read in to the neighbor
1050 write (iout,*) "Processor",me," entered WRITE_AND_SEND_CONF"
1053 call MPI_Send(icount,1,MPI_INTEGER,Next,570,MPI_COMM_WORLD,IERROR)
1054 call MPI_Send(nss_all(1),icount,MPI_INTEGER,&
1055 Next,571,MPI_COMM_WORLD,IERROR)
1056 call MPI_Send(ihpb_all(1,1),icount,MPI_INTEGER,&
1057 Next,572,MPI_COMM_WORLD,IERROR)
1058 call MPI_Send(jhpb_all(1,1),icount,MPI_INTEGER,&
1059 Next,573,MPI_COMM_WORLD,IERROR)
1060 call MPI_Send(rmstb(jj_old),icount,MPI_DOUBLE_PRECISION,&
1061 Next,577,MPI_COMM_WORLD,IERROR)
1062 call MPI_Send(entfac(jj_old),icount,MPI_DOUBLE_PRECISION,&
1063 Next,579,MPI_COMM_WORLD,IERROR)
1064 call MPI_Send(allcart(1,1,1),3*icount*2*nres,&
1065 MPI_REAL,Next,580,MPI_COMM_WORLD,IERROR)
1067 call dawrite_ccoords(jj_old,jj,icbase)
1069 end subroutine write_and_send_cconf
1070 !------------------------------------------------------------------------------
1072 subroutine receive_and_pass_cconf(icount,jj_old,jj,Previous,Next)
1076 ! include "DIMENSIONS"
1077 ! include "sizesclu.dat"
1079 integer :: IERROR !,STATUS(MPI_STATUS_SIZE)
1080 ! include "COMMON.MPI"
1081 ! include "COMMON.CHAIN"
1082 ! include "COMMON.SBRIDGE"
1083 ! include "COMMON.INTERACT"
1084 ! include "COMMON.IOUNITS"
1085 ! include "COMMON.VAR"
1086 ! include "COMMON.GEO"
1087 ! include "COMMON.CLUSTER"
1088 integer :: i,j,k,icount,jj_old,jj,Previous,Next
1091 write (iout,*) "Processor",me," entered RECEIVE_AND_PASS_CONF"
1094 do while (icount.gt.0)
1095 call MPI_Recv(icount,1,MPI_INTEGER,Previous,570,MPI_COMM_WORLD,&
1097 call MPI_Send(icount,1,MPI_INTEGER,Next,570,MPI_COMM_WORLD,&
1100 write (iout,*) "Processor",me," icount",icount
1102 if (icount.eq.0) return
1103 call MPI_Recv(nss_all(1),icount,MPI_INTEGER,&
1104 Previous,571,MPI_COMM_WORLD,STATUS,IERROR)
1105 call MPI_Send(nss_all(1),icount,MPI_INTEGER,&
1106 Next,571,MPI_COMM_WORLD,IERROR)
1107 call MPI_Recv(ihpb_all(1,1),icount,MPI_INTEGER,&
1108 Previous,572,MPI_COMM_WORLD,STATUS,IERROR)
1109 call MPI_Send(ihpb_all(1,1),icount,MPI_INTEGER,&
1110 Next,572,MPI_COMM_WORLD,IERROR)
1111 call MPI_Recv(jhpb_all(1,1),icount,MPI_INTEGER,&
1112 Previous,573,MPI_COMM_WORLD,STATUS,IERROR)
1113 call MPI_Send(jhpb_all(1,1),icount,MPI_INTEGER,&
1114 Next,573,MPI_COMM_WORLD,IERROR)
1115 call MPI_Recv(rmstb(jj_old),icount,MPI_DOUBLE_PRECISION,&
1116 Previous,577,MPI_COMM_WORLD,STATUS,IERROR)
1117 call MPI_Send(rmstb(jj_old),icount,MPI_DOUBLE_PRECISION,&
1118 Next,577,MPI_COMM_WORLD,IERROR)
1119 call MPI_Recv(entfac(jj_old),icount,MPI_DOUBLE_PRECISION,&
1120 Previous,579,MPI_COMM_WORLD,STATUS,IERROR)
1121 call MPI_Send(entfac(jj_old),icount,MPI_DOUBLE_PRECISION,&
1122 Next,579,MPI_COMM_WORLD,IERROR)
1123 call MPI_Recv(allcart(1,1,1),3*icount*2*nres,&
1124 MPI_REAL,Previous,580,MPI_COMM_WORLD,STATUS,IERROR)
1125 call MPI_Send(allcart(1,1,1),3*icount*2*nres,&
1126 MPI_REAL,Next,580,MPI_COMM_WORLD,IERROR)
1128 call dawrite_ccoords(jj_old,jj,icbase)
1131 write (iout,*) "Processor",me," received",icount," conformations"
1133 write (iout,'(8f10.4)') (allcart(l,k,i),l=1,3,k=1,nres)
1134 write (iout,'(8f10.4)')((allcart(l,k,i+nres),l=1,3,k=nnt,nct)
1135 write (iout,'(e15.5,16i5)') entfac(i)
1140 end subroutine receive_and_pass_cconf
1142 !------------------------------------------------------------------------------
1143 subroutine daread_ccoords(istart_conf,iend_conf)
1145 use energy_data, only: dyn_ss
1148 ! include "DIMENSIONS"
1149 ! include "sizesclu.dat"
1153 ! include "COMMON.MPI"
1155 ! include "COMMON.CHAIN"
1156 ! include "COMMON.CLUSTER"
1157 ! include "COMMON.IOUNITS"
1158 ! include "COMMON.INTERACT"
1159 ! include "COMMON.VAR"
1160 ! include "COMMON.SBRIDGE"
1161 ! include "COMMON.GEO"
1162 integer :: istart_conf,iend_conf
1163 integer :: i,j,ij,ii,iii
1165 character(len=16) :: form,acc
1166 character(len=32) :: nam
1168 ! Read conformations off a DA scratchfile.
1171 write (iout,*) "DAREAD_COORDS"
1172 write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
1173 inquire(unit=icbase,name=nam,recl=len,form=form,access=acc)
1174 write (iout,*) "len=",len," form=",form," acc=",acc
1175 write (iout,*) "nam=",nam
1178 do ii=istart_conf,iend_conf
1179 ij = ii - istart_conf + 1
1182 write (iout,*) "Reading binary file, record",iii," ii",ii
1186 read(icbase,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres), &
1187 ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
1188 entfac(ii),rmstb(ii)
1190 read(icbase,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres),&
1191 ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
1192 nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss),&
1193 entfac(ii),rmstb(ii)
1197 write (iout,*) ii,iii,ij,entfac(ii)
1198 write (iout,'(8f10.5)') ((allcart(j,i,ij),j=1,3),i=1,nres)
1199 write (iout,'(8f10.4)') ((allcart(j,i,ij),j=1,3),&
1200 i=nnt+nres,nct+nres)
1201 write (iout,'(2e15.5)') entfac(ij)
1202 write (iout,'(16i5)') nss_all(ij),(ihpb_all(i,ij),&
1203 jhpb_all(i,ij),i=1,nss)
1208 end subroutine daread_ccoords
1209 !------------------------------------------------------------------------------
1210 subroutine dawrite_ccoords(istart_conf,iend_conf,unit_out)
1213 ! include "DIMENSIONS"
1214 ! include "sizesclu.dat"
1215 use energy_data, only: dyn_ss
1220 ! include "COMMON.MPI"
1222 ! include "COMMON.CHAIN"
1223 ! include "COMMON.INTERACT"
1224 ! include "COMMON.IOUNITS"
1225 ! include "COMMON.VAR"
1226 ! include "COMMON.SBRIDGE"
1227 ! include "COMMON.GEO"
1228 ! include "COMMON.CLUSTER"
1229 integer :: istart_conf,iend_conf
1230 integer :: i,j,ii,ij,iii,unit_out
1232 character(len=16) :: form,acc
1233 character(len=32) :: nam
1235 ! Write conformations to a DA scratchfile.
1238 write (iout,*) "DAWRITE_COORDS"
1239 write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
1240 write (iout,*) "lenrec",lenrec
1241 inquire(unit=unit_out,name=nam,recl=len,form=form,access=acc)
1242 write (iout,*) "len=",len," form=",form," acc=",acc
1243 write (iout,*) "nam=",nam
1246 do ii=istart_conf,iend_conf
1248 ij = ii - istart_conf + 1
1250 write (iout,*) "Writing binary file, record",iii," ii",ii
1254 write(unit_out,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres),&
1255 ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
1256 entfac(ii),rmstb(ii)
1258 write(unit_out,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres), &
1259 ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
1260 nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss_all(ij)),&
1261 entfac(ii),rmstb(ii)
1265 write (iout,'(8f10.5)') ((allcart(j,i,ij),j=1,3),i=1,nres)
1266 write (iout,'(8f10.4)') ((allcart(j,i,ij),j=1,3),i=nnt+nres,&
1268 write (iout,'(2e15.5)') entfac(ij)
1269 write (iout,'(16i5)') nss_all(ij),(ihpb(i,ij),jhpb(i,ij),i=1,&
1275 end subroutine dawrite_ccoords
1276 !-----------------------------------------------------------------------------
1278 !-----------------------------------------------------------------------------
1279 subroutine read_control
1281 ! Read molecular data
1283 use energy_data, only: rescale_mode,distchainmax,ipot,constr_homology,&
1284 with_dihed_constr,read_homol_frag,out_template_coord,&
1286 use control_data, only: titel,outpdb,outmol2,refstr,pdbref,&
1287 iscode,symetr,punch_dist,print_dist,nstart,nend,&
1288 caonly,iopt,efree,lprint_cart,lprint_int,rlamb_ele,&
1289 r_cut_ele,nclust,tor_mode,scelemode,r_cut_mart,r_cut_ang,&
1290 rlamb_mart,constr_dist
1291 use geometry_data, only:bordliptop,bordlipbot,&
1292 bufliptop,buflipbot,lipthick,lipbufthick
1294 ! include 'DIMENSIONS'
1295 ! include 'sizesclu.dat'
1296 ! include 'COMMON.IOUNITS'
1297 ! include 'COMMON.TIME1'
1298 ! include 'COMMON.SBRIDGE'
1299 ! include 'COMMON.CONTROL'
1300 ! include 'COMMON.CLUSTER'
1301 ! include 'COMMON.CHAIN'
1302 ! include 'COMMON.HEADER'
1303 ! include 'COMMON.FFIELD'
1304 ! include 'COMMON.FREE'
1305 ! include 'COMMON.INTERACT'
1306 character(len=320) :: controlcard !,ucase
1308 ! include 'COMMON.INFO'
1312 read (INP,'(a80)') titel
1313 call card_concat(controlcard,.true.)
1315 call readi(controlcard,'NRES',nres_molec(1),0)
1316 call readi(controlcard,'NRES_NUCL',nres_molec(2),0)
1317 call readi(controlcard,'NRES_CAT',nres_molec(5),0)
1320 nres=nres_molec(i)+nres
1322 print *,"TU",nres_molec(:)
1324 ! call alloc_clust_arrays
1325 allocate(rcutoff(max_cut+1)) !(max_cut+1)
1326 allocate(beta_h(maxT)) !(maxT)
1327 allocate(mult(nres)) !(maxres)
1330 call readi(controlcard,'RESCALE',rescale_mode,2)
1331 call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
1332 write (iout,*) "DISTCHAINMAX",distchainmax
1333 call readi(controlcard,'PDBOUT',outpdb,0)
1334 call readi(controlcard,'MOL2OUT',outmol2,0)
1335 refstr=(index(controlcard,'REFSTR').gt.0)
1336 write (iout,*) "REFSTR",refstr
1337 pdbref=(index(controlcard,'PDBREF').gt.0)
1338 iscode=index(controlcard,'ONE_LETTER')
1339 tree=(index(controlcard,'MAKE_TREE').gt.0)
1340 min_var=(index(controlcard,'MINVAR').gt.0)
1341 plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
1342 punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
1343 call readi(controlcard,'NCUT',ncut,1)
1344 call readi(controlcard,'SYM',symetr,1)
1345 write (iout,*) 'sym', symetr
1346 call readi(controlcard,'NSTART',nstart,0)
1347 call reada(controlcard,'BOXX',boxxsize,100.0d0)
1348 call reada(controlcard,'BOXY',boxysize,100.0d0)
1349 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
1350 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
1351 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
1352 if (lipthick.gt.0.0d0) then
1353 bordliptop=(boxzsize+lipthick)/2.0
1354 bordlipbot=bordliptop-lipthick
1355 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
1356 write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
1357 buflipbot=bordlipbot+lipbufthick
1358 bufliptop=bordliptop-lipbufthick
1359 if ((lipbufthick*2.0d0).gt.lipthick) &
1360 write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
1361 endif !lipthick.gt.0
1362 write(iout,*) "bordliptop=",bordliptop
1363 write(iout,*) "bordlipbot=",bordlipbot
1364 write(iout,*) "bufliptop=",bufliptop
1365 write(iout,*) "buflipbot=",buflipbot
1367 call readi(controlcard,'NCLUST',nclust,5)
1368 ! ions=index(controlcard,"IONS").gt.0
1369 call readi(controlcard,"SCELEMODE",scelemode,0)
1370 print *,"SCELE",scelemode
1371 call readi(controlcard,'TORMODE',tor_mode,0)
1372 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1373 write(iout,*) "torsional and valence angle mode",tor_mode
1374 call reada(controlcard,"R_CUT_ELE",r_cut_ele,25.0d0)
1375 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
1376 write(iout,*) "R_CUT_ELE=",r_cut_ele
1377 call reada(controlcard,"R_CUT_MART",r_cut_mart,15.0d0)
1378 call reada(controlcard,"LAMBDA_MART",rlamb_mart,0.3d0)
1379 call reada(controlcard,"R_CUT_ANG",r_cut_ang,4.2d0)
1380 call readi(controlcard,'NEND',nend,0)
1381 call reada(controlcard,'ECUT',ecut,10.0d0)
1382 call reada(controlcard,'PROB',prob_limit,0.99d0)
1383 write (iout,*) "Probability limit",prob_limit
1384 lgrp=(index(controlcard,'LGRP').gt.0)
1385 caonly=(index(controlcard,'CA_ONLY').gt.0)
1386 print_dist=(index(controlcard,'PRINT_DIST').gt.0)
1387 call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
1388 call readi(controlcard,'IOPT',iopt,2)
1389 lside = index(controlcard,"SIDE").gt.0
1390 efree = index(controlcard,"EFREE").gt.0
1391 call readi(controlcard,'NTEMP',nT,1)
1392 write (iout,*) "nT",nT
1393 !el call reada(controlcard,'TEMP0',temp0,300.0d0) !el
1394 call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
1395 write (iout,*) "nT",nT
1396 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
1398 beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
1400 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
1401 lprint_cart=index(controlcard,"PRINT_CART") .gt.0
1402 lprint_int=index(controlcard,"PRINT_INT") .gt.0
1403 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
1404 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
1405 !c if (constr_homology) tole=dmax1(tole,1.5d0)
1406 write (iout,*) "with_homology_constr ",with_dihed_constr,&
1407 " CONSTR_HOMOLOGY",constr_homology
1408 read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
1409 out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0
1410 out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0
1411 write (iout,*) "out_template_coord ",OUT_TEMPLATE_COORD
1414 end subroutine read_control
1415 !-----------------------------------------------------------------------------
1416 subroutine read_constr_homology
1421 use control, only:init_int_table,homology_partition
1422 use MD_data, only:iset
1424 ! include 'DIMENSIONS'
1428 ! include 'COMMON.SETUP'
1429 ! include 'COMMON.CONTROL'
1430 ! include 'COMMON.HOMOLOGY'
1431 ! include 'COMMON.CHAIN'
1432 ! include 'COMMON.IOUNITS'
1433 ! include 'COMMON.MD'
1434 ! include 'COMMON.QRESTR'
1435 ! include 'COMMON.GEO'
1436 ! include 'COMMON.INTERACT'
1437 ! include 'COMMON.NAMES'
1438 ! include 'COMMON.VAR'
1441 ! double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
1443 ! common /przechowalnia/ odl_temp(maxres,maxres,max_template),
1444 ! & sigma_odl_temp(maxres,maxres,max_template)
1446 character*24 model_ki_dist, model_ki_angle
1447 character*500 controlcard
1448 integer :: ki,i,ii,j,k,l
1449 integer, dimension (:), allocatable :: ii_in_use
1450 integer :: i_tmp,idomain_tmp,&
1451 irec,ik,iistart,nres_temp
1454 logical :: liiflag,lfirst
1457 ! FP - Nov. 2014 Temporary specifications for new vars
1459 real(kind=8) :: rescore_tmp,x12,y12,z12,rescore2_tmp,&
1460 rescore3_tmp, dist_cut
1461 real(kind=8), dimension (:,:),allocatable :: rescore
1462 real(kind=8), dimension (:,:),allocatable :: rescore2
1463 real(kind=8), dimension (:,:),allocatable :: rescore3
1464 real(kind=8) :: distal
1465 character*24 tpl_k_rescore
1466 character*256 pdbfile
1468 ! -----------------------------------------------------------------
1469 ! Reading multiple PDB ref structures and calculation of retraints
1470 ! not using pre-computed ones stored in files model_ki_{dist,angle}
1472 ! -----------------------------------------------------------------
1475 ! Alternative: reading from input
1476 call card_concat(controlcard,.true.)
1477 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
1478 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
1479 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
1480 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
1481 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
1482 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
1483 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
1484 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
1485 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
1486 if(.not.read2sigma.and.start_from_model) then
1487 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)&
1488 write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
1489 start_from_model=.false.
1490 iranconf=(indpdb.le.0)
1492 if(start_from_model .and. (me.eq.king .or. .not. out1file))&
1493 write(iout,*) 'START_FROM_MODELS is ON'
1494 ! if(start_from_model .and. rest) then
1495 ! if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1496 ! write(iout,*) 'START_FROM_MODELS is OFF'
1497 ! write(iout,*) 'remove restart keyword from input'
1500 if (start_from_model) nmodel_start=constr_homology
1501 if(.not.allocated(waga_homology)) allocate (waga_homology(homol_nset))
1502 if (homol_nset.gt.1)then
1503 call card_concat(controlcard,.true.)
1504 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
1505 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1506 ! write(iout,*) "iset homology_weight "
1508 write(iout,*) i,waga_homology(i)
1511 iset=mod(kolor,homol_nset)+1
1514 waga_homology(1)=1.0
1517 !d write (iout,*) "nnt",nnt," nct",nct
1520 if (read_homol_frag) then
1521 call read_klapaucjusz
1527 ! write(iout,*) 'nnt=',nnt,'nct=',nct
1530 ! do k=1,constr_homology
1544 if(.not.allocated(pdbfiles_chomo)) allocate(pdbfiles_chomo(constr_homology))
1545 if(.not.allocated(chomo)) allocate(chomo(3,2*nres+2,constr_homology))
1547 do k=1,constr_homology
1549 read(inp,'(a)') pdbfile
1550 pdbfiles_chomo(k)=pdbfile
1551 if(me.eq.king .or. .not. out1file) &
1552 write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',&
1553 pdbfile(:ilen(pdbfile))
1554 open(ipdbin,file=pdbfile,status='old',err=33)
1556 33 write (iout,'(a,5x,a)') 'Error opening PDB file',&
1557 pdbfile(:ilen(pdbfile))
1560 ! print *,'Begin reading pdb data'
1562 ! Files containing res sim or local scores (former containing sigmas)
1565 write(kic2,'(bz,i2.2)') k
1567 tpl_k_rescore="template"//kic2//".sco"
1568 write(iout,*) "tpl_k_rescore=",tpl_k_rescore
1571 write(iout,*) "read2sigma",read2sigma
1573 if (read2sigma) then
1574 call readpdb_template(k)
1578 write(iout,*) "after readpdb"
1579 if(.not.allocated(nres_chomo)) allocate(nres_chomo(constr_homology))
1582 if(.not.allocated(rescore)) allocate(rescore(constr_homology,nres))
1583 if(.not.allocated(rescore2)) allocate(rescore2(constr_homology,nres))
1584 if(.not.allocated(rescore3)) allocate(rescore3(constr_homology,nres))
1585 if(.not.allocated(ii_in_use)) allocate(ii_in_use(nres*(nres-1)))
1586 if(.not.allocated(idomain)) allocate(idomain(constr_homology,nres))
1587 if(.not.allocated(l_homo)) allocate(l_homo(constr_homology,1000*nres))
1588 if(.not.allocated(ires_homo)) allocate(ires_homo(nres*200))
1589 if(.not.allocated(jres_homo)) allocate(jres_homo(nres*200))
1590 if(.not.allocated(odl)) allocate(odl(constr_homology,nres*200))
1591 if(.not.allocated(sigma_odl)) allocate(sigma_odl(constr_homology,nres*200))
1592 if(.not.allocated(dih)) allocate(dih(constr_homology,nres))
1593 if(.not.allocated(sigma_dih)) allocate(sigma_dih(constr_homology,nres))
1594 if(.not.allocated(thetatpl)) allocate(thetatpl(constr_homology,nres))
1595 if(.not.allocated(sigma_theta)) allocate(sigma_theta(constr_homology,nres))
1596 ! if(.not.allocated(thetatpl)) allocate(thetatpl(constr_homology,nres))
1597 if(.not.allocated(sigma_d)) allocate(sigma_d(constr_homology,nres))
1598 if(.not.allocated(xxtpl)) allocate(xxtpl(constr_homology,nres))
1599 if(.not.allocated(yytpl)) allocate(yytpl(constr_homology,nres))
1600 if(.not.allocated(zztpl)) allocate(zztpl(constr_homology,nres))
1601 ! if(.not.allocated(distance)) allocate(distance(constr_homology))
1602 ! if(.not.allocated(distancek)) allocate(distancek(constr_homology))
1606 ! Distance restraints
1609 ! Copy the coordinates from reference coordinates (?)
1610 do i=1,2*nres_chomo(k)
1613 ! write (iout,*) "c(",j,i,") =",c(j,i)
1617 ! From read_dist_constr (commented out 25/11/2014 <-> res sim)
1619 ! write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
1620 open (ientin,file=tpl_k_rescore,status='old')
1621 if (nnt.gt.1) rescore(k,1)=0.0d0
1622 do irec=nnt,nct ! loop for reading res sim
1623 if (read2sigma) then
1624 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,&
1625 rescore3_tmp,idomain_tmp
1627 idomain(k,i_tmp)=idomain_tmp
1628 rescore(k,i_tmp)=rescore_tmp
1629 rescore2(k,i_tmp)=rescore2_tmp
1630 rescore3(k,i_tmp)=rescore3_tmp
1631 if (.not. out1file .or. me.eq.king)&
1632 write(iout,'(a7,i5,3f10.5,i5)') "rescore",&
1633 i_tmp,rescore2_tmp,rescore_tmp,&
1634 rescore3_tmp,idomain_tmp
1637 read (ientin,*,end=1401) rescore_tmp
1639 ! rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
1640 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
1641 ! write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
1646 if (waga_dist.ne.0.0d0) then
1654 distal=dsqrt(x12*x12+y12*y12+z12*z12)
1655 ! write (iout,*) k,i,j,distal,dist2_cut
1657 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 &
1658 .and. distal.le.dist2_cut ) then
1664 ! write (iout,*) "k",k
1665 ! write (iout,*) "i",i," j",j," constr_homology",
1670 if (read2sigma) then
1673 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
1675 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
1676 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = &
1677 sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1679 if (odl(k,ii).le.dist_cut) then
1680 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
1683 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
1684 dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
1686 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
1687 dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
1691 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
1694 ! l_homo(k,ii)=.false.
1700 ! write (iout,*) "Distance restraints set"
1703 ! Theta, dihedral and SC retraints
1705 if (waga_angle.gt.0.0d0) then
1706 ! open (ientin,file=tpl_k_sigma_dih,status='old')
1707 ! do irec=1,maxres-3 ! loop for reading sigma_dih
1708 ! read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
1709 ! if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
1710 ! sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
1711 ! & sigma_dih(k,i+nnt-1)
1716 if (idomain(k,i).eq.0) then
1720 dih(k,i)=phiref(i) ! right?
1721 ! read (ientin,*) sigma_dih(k,i) ! original variant
1722 ! write (iout,*) "dih(",k,i,") =",dih(k,i)
1723 ! write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1724 ! & "rescore(",k,i-1,") =",rescore(k,i-1),
1725 ! & "rescore(",k,i-2,") =",rescore(k,i-2),
1726 ! & "rescore(",k,i-3,") =",rescore(k,i-3)
1728 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+ &
1729 rescore(k,i-2)+rescore(k,i-3))/4.0
1730 ! if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
1731 ! write (iout,*) "Raw sigmas for dihedral angle restraints"
1732 ! write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
1733 ! sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1734 ! rescore(k,i-2)*rescore(k,i-3) ! right expression ?
1735 ! Instead of res sim other local measure of b/b str reliability possible
1736 if (sigma_dih(k,i).ne.0) &
1737 sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
1738 ! sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
1742 ! write (iout,*) "Dihedral angle restraints set"
1745 if (waga_theta.gt.0.0d0) then
1746 ! open (ientin,file=tpl_k_sigma_theta,status='old')
1747 ! do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
1748 ! read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
1749 ! sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
1750 ! & sigma_theta(k,i+nnt-1)
1755 do i = nnt+2,nct ! right? without parallel.
1756 ! do i = i=1,nres ! alternative for bounds acc to readpdb?
1757 ! do i=ithet_start,ithet_end ! with FG parallel.
1758 if (idomain(k,i).eq.0) then
1759 sigma_theta(k,i)=0.0
1762 thetatpl(k,i)=thetaref(i)
1763 ! write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
1764 ! write(iout,*) "rescore(",k,i,") =",rescore(k,i),
1765 ! & "rescore(",k,i-1,") =",rescore(k,i-1),
1766 ! & "rescore(",k,i-2,") =",rescore(k,i-2)
1767 ! read (ientin,*) sigma_theta(k,i) ! 1st variant
1768 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+ &
1770 ! if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
1771 if (sigma_theta(k,i).ne.0) &
1772 sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
1774 ! sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
1775 ! rescore(k,i-2) ! right expression ?
1776 ! sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
1779 ! write (iout,*) "Angle restraints set"
1782 if (waga_d.gt.0.0d0) then
1783 ! open (ientin,file=tpl_k_sigma_d,status='old')
1784 ! do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
1785 ! read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
1786 ! sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
1787 ! & sigma_d(k,i+nnt-1)
1791 do i = nnt,nct ! right? without parallel.
1792 ! do i=2,nres-1 ! alternative for bounds acc to readpdb?
1793 ! do i=loc_start,loc_end ! with FG parallel.
1794 if (itype(i,1).eq.10) cycle
1795 if (idomain(k,i).eq.0 ) then
1802 ! write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
1803 ! write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
1804 ! write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
1805 ! write(iout,*) "rescore(",k,i,") =",rescore(k,i)
1806 sigma_d(k,i)=rescore3(k,i) ! right expression ?
1807 if (sigma_d(k,i).ne.0) &
1808 sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
1810 ! sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
1811 ! sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
1812 ! read (ientin,*) sigma_d(k,i) ! 1st variant
1816 ! write (iout,*) "SC restraints set"
1819 ! remove distance restraints not used in any model from the list
1820 ! shift data in all arrays
1822 ! write (iout,*) "waga_dist",waga_dist," nnt",nnt," nct",nct
1823 if (waga_dist.ne.0.0d0) then
1830 ! if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
1831 ! & .and. distal.le.dist2_cut ) then
1832 ! write (iout,*) "i",i," j",j," ii",ii
1834 if (ii_in_use(ii).eq.0.and.liiflag.or. &
1835 ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
1842 if(i10.eq.lim_odl) i10=i10+1
1844 ires_homo(iistart+ki)=ires_homo(ki+i01)
1845 jres_homo(iistart+ki)=jres_homo(ki+i01)
1846 ii_in_use(iistart+ki)=ii_in_use(ki+i01)
1847 do k=1,constr_homology
1848 odl(k,iistart+ki)=odl(k,ki+i01)
1849 sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
1850 l_homo(k,iistart+ki)=l_homo(k,ki+i01)
1853 iistart=iistart+i10-i01
1856 if (ii_in_use(ii).ne.0.and..not.liiflag) then
1864 ! write (iout,*) "Removing distances completed"
1866 endif ! .not. klapaucjusz
1868 if (constr_homology.gt.0) call homology_partition
1869 write (iout,*) "After homology_partition"
1871 if (constr_homology.gt.0) call init_int_table
1872 write (iout,*) "After init_int_table"
1874 ! endif ! .not. klapaucjusz
1876 ! if (constr_homology.gt.0) call homology_partition
1877 ! write (iout,*) "After homology_partition"
1879 ! if (constr_homology.gt.0) call init_int_table
1880 ! write (iout,*) "After init_int_table"
1882 ! write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
1883 ! write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
1887 if (.not.out_template_restr) return
1888 !d write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
1889 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1890 write (iout,*) "Distance restraints from templates"
1892 write(iout,'(3i7,100(2f8.2,1x,l1,4x))') &
1893 ii,ires_homo(ii),jres_homo(ii),&
1894 (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),&
1895 ki=1,constr_homology)
1897 write (iout,*) "Dihedral angle restraints from templates"
1899 write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i,1),1),&
1900 (rad2deg*dih(ki,i),&
1901 rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
1903 write (iout,*) "Virtual-bond angle restraints from templates"
1905 write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i,1),1),&
1906 (rad2deg*thetatpl(ki,i),&
1907 rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
1909 write (iout,*) "SC restraints from templates"
1911 write(iout,'(i7,100(4f8.2,4x))') i,&
1912 (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i), &
1913 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
1917 end subroutine read_constr_homology
1918 !----------------------------------------------------------------------------
1921 ! Read molecular data.
1923 use geometry_data, only: nsup,cref,nres0,nstart_sup,nstart_seq,dc,&
1924 deg2rad,phibound,crefjlee,cref,rad2deg
1925 use energy_data!, only: wsc,wscp,welec,wbond,wstrain,wtor,wtor_d,&
1926 ! wang,wscloc,wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,&
1927 ! wturn3,wturn4,wturn6,wvdwpp,weights
1928 use control_data, only: titel,nstart,nend,pdbref,refstr,iscode,&
1929 indpdb,constr_dist,raw_psipred, with_theta_constr
1930 use geometry, only: chainbuild,alloc_geo_arrays
1931 use energy, only: alloc_ener_arrays
1932 use io_base, only: rescode
1933 use control, only: setup_var,init_int_table
1934 use conform_compar, only: contact
1936 ! include 'DIMENSIONS'
1937 ! include 'COMMON.IOUNITS'
1938 ! include 'COMMON.GEO'
1939 ! include 'COMMON.VAR'
1940 ! include 'COMMON.INTERACT'
1941 ! include 'COMMON.LOCAL'
1942 ! include 'COMMON.NAMES'
1943 ! include 'COMMON.CHAIN'
1944 ! include 'COMMON.FFIELD'
1945 ! include 'COMMON.SBRIDGE'
1946 ! include 'COMMON.HEADER'
1947 ! include 'COMMON.CONTROL'
1948 ! include 'COMMON.CONTACTS'
1949 ! include 'COMMON.TIME1'
1951 ! include 'COMMON.INFO'
1953 character(len=4) :: sequence(nres) !(maxres)
1954 character(len=800) :: weightcard,controlcard
1956 real(kind=8) :: x(6*nres) !(maxvar)
1957 real(kind=8) :: ssscale
1958 real(kind=8) :: phihel,phibet,sigmahel,sigmabet,sumv,&
1960 integer :: itype_pdb(nres) !(maxres)
1962 integer :: i,j,kkk,mnum
1966 !el allocate(weights(n_ene))
1967 allocate(weights(max_ene))
1968 call alloc_geo_arrays
1969 call alloc_ener_arrays
1970 !-----------------------------
1971 allocate(c(3,2*nres+2)) !(3,maxres2+2) maxres2=2*maxres
1972 allocate(dc(3,0:2*nres+2)) !(3,0:maxres2)
1973 allocate(itype(nres+2,5)) !(maxres)
1974 allocate(molnum(nres+1))
1975 allocate(itel(nres+2))
1989 !--------------------------
1990 ! Read weights of the subsequent energy terms.
1991 call card_concat(weightcard,.true.)
1992 call reada(weightcard,'WSC',wsc,1.0d0)
1993 call reada(weightcard,'WLONG',wsc,wsc)
1994 call reada(weightcard,'WSCP',wscp,1.0d0)
1995 call reada(weightcard,'WELEC',welec,1.0D0)
1996 call reada(weightcard,'WVDWPP',wvdwpp,welec)
1997 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
1998 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
1999 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
2000 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
2001 call reada(weightcard,'WTURN3',wturn3,1.0D0)
2002 call reada(weightcard,'WTURN4',wturn4,1.0D0)
2003 call reada(weightcard,'WTURN6',wturn6,1.0D0)
2004 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
2005 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
2006 call reada(weightcard,'WBOND',wbond,1.0D0)
2007 call reada(weightcard,'WTOR',wtor,1.0D0)
2008 call reada(weightcard,'WTORD',wtor_d,1.0D0)
2009 call reada(weightcard,'WANG',wang,1.0D0)
2010 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
2011 call reada(weightcard,'SCAL14',scal14,0.4D0)
2012 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
2013 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
2014 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
2015 call reada(weightcard,'WCATCAT',wcatcat,0.0d0)
2016 call reada(weightcard,'WCATPROT',wcatprot,0.0d0)
2017 call reada(weightcard,'TEMP0',temp0,300.0d0) !!! el
2018 if (index(weightcard,'SOFT').gt.0) ipot=6
2019 call reada(weightcard,'WVDWPP_NUCL',wvdwpp_nucl,0.0D0)
2020 call reada(weightcard,'WELPP',welpp,0.0d0)
2021 call reada(weightcard,'WVDWPSB',wvdwpsb,0.0d0)
2022 call reada(weightcard,'WELPSB',welpsb,0.0D0)
2023 call reada(weightcard,'WVDWSB',wvdwsb,0.0d0)
2024 call reada(weightcard,'WELSB',welsb,0.0D0)
2025 call reada(weightcard,'WBOND_NUCL',wbond_nucl,0.0D0)
2026 call reada(weightcard,'WANG_NUCL',wang_nucl,0.0D0)
2027 call reada(weightcard,'WSBLOC',wsbloc,0.0D0)
2028 call reada(weightcard,'WTOR_NUCL',wtor_nucl,0.0D0)
2029 ! print *,"KUR..",wtor_nucl
2030 call reada(weightcard,'WTORD_NUCL',wtor_d_nucl,0.0D0)
2031 call reada(weightcard,'WCORR_NUCL',wcorr_nucl,0.0D0)
2032 call reada(weightcard,'WCORR3_NUC',wcorr3_nucl,0.0D0)
2033 ! 12/1/95 Added weight for the multi-body term WCORR
2034 call reada(weightcard,'WCORRH',wcorr,1.0D0)
2035 call reada(weightcard,'WSCBASE',wscbase,0.0D0)
2036 call reada(weightcard,'WPEPBASE',wpepbase,1.0d0)
2037 call reada(weightcard,'WSCPHO',wscpho,0.0d0)
2038 call reada(weightcard,'WPEPPHO',wpeppho,0.0d0)
2039 call reada(weightcard,'WCATNUCL',wcatnucl,0.0d0)
2040 call reada(weightcard,'WCATTRAN',wcat_tran,0.0d0)
2041 call reada(weightcard,"D0CM",d0cm,3.78d0)
2042 call reada(weightcard,"AKCM",akcm,15.1d0)
2043 call reada(weightcard,"AKTH",akth,11.0d0)
2044 call reada(weightcard,"AKCT",akct,12.0d0)
2045 call reada(weightcard,"V1SS",v1ss,-1.08d0)
2046 call reada(weightcard,"V2SS",v2ss,7.61d0)
2047 call reada(weightcard,"V3SS",v3ss,13.7d0)
2048 call reada(weightcard,"EBR",ebr,-5.50D0)
2049 call reada(weightcard,"ATRISS",atriss,0.301D0)
2050 call reada(weightcard,"BTRISS",btriss,0.021D0)
2051 call reada(weightcard,"CTRISS",ctriss,1.001D0)
2052 call reada(weightcard,"DTRISS",dtriss,1.001D0)
2053 call reada(weightcard,"SSSCALE",ssscale,1.0D0)
2054 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
2055 call reada(weightcard,"LIPSCALE",lipscale,1.0D0)
2056 call reada(weightcard,"WTL",wliptran,1.0D0)
2058 call reada(weightcard,"HT",Ht,0.0D0)
2060 ss_depth=(ebr/wsc-0.25*eps(1,1))*ssscale
2061 Ht=(Ht/wsc-0.25*eps(1,1))*ssscale
2062 akcm=akcm/wsc*ssscale
2063 akth=akth/wsc*ssscale
2064 akct=akct/wsc*ssscale
2065 v1ss=v1ss/wsc*ssscale
2066 v2ss=v2ss/wsc*ssscale
2067 v3ss=v3ss/wsc*ssscale
2069 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
2072 if (wcorr4.gt.0.0d0) wcorr=wcorr4
2091 !el weights(19)=wsccor !!!!!!!!!!!!!!!!
2093 weights(26)=wvdwpp_nucl
2099 weights(32)=wbond_nucl
2100 weights(33)=wang_nucl
2102 weights(35)=wtor_nucl
2103 weights(36)=wtor_d_nucl
2104 weights(37)=wcorr_nucl
2105 weights(38)=wcorr3_nucl
2107 weights(42)=wcatprot
2112 weights(50)=wcatnucl
2113 weights(56)=wcat_tran
2114 weights(22)=wliptran
2116 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
2117 wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,&
2118 wturn4,wturn6,wsccor
2119 10 format (/'Energy-term weights (unscaled):'// &
2120 'WSCC= ',f10.6,' (SC-SC)'/ &
2121 'WSCP= ',f10.6,' (SC-p)'/ &
2122 'WELEC= ',f10.6,' (p-p electr)'/ &
2123 'WVDWPP= ',f10.6,' (p-p VDW)'/ &
2124 'WBOND= ',f10.6,' (stretching)'/ &
2125 'WANG= ',f10.6,' (bending)'/ &
2126 'WSCLOC= ',f10.6,' (SC local)'/ &
2127 'WTOR= ',f10.6,' (torsional)'/ &
2128 'WTORD= ',f10.6,' (double torsional)'/ &
2129 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ &
2130 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ &
2131 'WCORR4= ',f10.6,' (multi-body 4th order)'/ &
2132 'WCORR5= ',f10.6,' (multi-body 5th order)'/ &
2133 'WCORR6= ',f10.6,' (multi-body 6th order)'/ &
2134 'WTURN3= ',f10.6,' (turns, 3rd order)'/ &
2135 'WTURN4= ',f10.6,' (turns, 4th order)'/ &
2136 'WTURN6= ',f10.6,' (turns, 6th order)'/ &
2137 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
2139 if (wcorr4.gt.0.0d0) then
2140 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',&
2141 'between contact pairs of peptide groups'
2142 write (iout,'(2(a,f5.3/))') &
2143 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,&
2144 'Range of quenching the correlation terms:',2*delt_corr
2145 else if (wcorr.gt.0.0d0) then
2146 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',&
2147 'between contact pairs of peptide groups'
2149 write (iout,'(a,f8.3)') &
2150 'Scaling factor of 1,4 SC-p interactions:',scal14
2151 write (iout,'(a,f8.3)') &
2152 'General scaling factor of SC-p interactions:',scalscp
2153 r0_corr=cutoff_corr-delt_corr
2155 aad(i,1)=scalscp*aad(i,1)
2156 aad(i,2)=scalscp*aad(i,2)
2157 bad(i,1)=scalscp*bad(i,1)
2158 bad(i,2)=scalscp*bad(i,2)
2162 print *,'indpdb=',indpdb,' pdbref=',pdbref
2164 ! Read sequence if not taken from the pdb file.
2165 if (iscode.gt.0) then
2166 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
2168 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
2170 print *,nres_molec(:),nres
2171 ! Convert sequence to numeric code
2172 do i=1,nres_molec(1)
2175 itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
2177 do i=nres_molec(1)+1,nres_molec(1)+nres_molec(2)
2180 write (iout,*),i,sequence(i)
2181 itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
2184 do i=nres_molec(1)+nres_molec(2)+1,nres_molec(1)+nres_molec(2)+nres_molec(5)
2187 itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
2190 print '(20i4)',(itype(i,molnum(i)),i=1,nres)
2194 if (itype(i,molnum(i)).eq.ntyp1_molec(molnum(i))&
2195 .or. itype(i+1,molnum(i+1)).eq.ntyp1_molec(molnum(i+1))) then
2197 if (itype(i,molnum(i)).eq.ntyp1_molec(molnum(i))) then
2201 else if (iabs(itype(i+1,1)).ne.20) then
2203 else if (iabs(itype(i,1)).ne.20) then
2210 write (iout,*) "ITEL"
2212 write (iout,*) i,itype(i,molnum(i)),itel(i)
2215 print *,'Call Read_Bridge.'
2219 print *,'NNT=',NNT,' NCT=',NCT
2220 if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2
2221 if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1
2222 if (nstart.lt.nnt) nstart=nnt
2223 if (nend.gt.nct .or. nend.eq.0) nend=nct
2224 write (iout,*) "nstart",nstart," nend",nend
2227 ! read(inp,'(a)') pdbfile
2228 ! write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
2229 ! open(ipdbin,file=pdbfile,status='old',err=33)
2231 ! 33 write (iout,'(a)') 'Error opening PDB file.'
2234 ! print *,'Begin reading pdb data'
2236 ! print *,'Finished reading pdb data'
2237 ! write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
2239 ! itype_pdb(i)=itype(i)
2242 ! write (iout,'(a,i3)') 'nsup=',nsup
2244 ! if (nsup.le.(nct-nnt+1)) then
2245 ! do i=0,nct-nnt+1-nsup
2246 ! if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
2251 ! write (iout,'(a)')
2252 ! & 'Error - sequences to be superposed do not match.'
2255 ! do i=0,nsup-(nct-nnt+1)
2256 ! if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
2258 ! nstart_sup=nstart_sup+i
2263 ! write (iout,'(a)')
2264 ! & 'Error - sequences to be superposed do not match.'
2267 ! write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
2268 ! & ' nstart_seq=',nstart_seq
2270 if (with_dihed_constr) then
2272 read (inp,*) ndih_constr
2273 if (ndih_constr.gt.0) then
2275 !C read (inp,*) ftors
2276 !C write (iout,*) 'FTORS',ftors
2277 !C ftors is the force constant for torsional quartic constrains
2278 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),&
2281 'There are',ndih_constr,' constraints on phi angles.'
2283 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),&
2287 phi0(i)=deg2rad*phi0(i)
2288 drange(i)=deg2rad*drange(i)
2290 else if (ndih_constr.lt.0) then
2292 call card_concat(controlcard,.true.)
2293 call reada(controlcard,"PHIHEL",phihel,50.0D0)
2294 call reada(controlcard,"PHIBET",phibet,180.0D0)
2295 call reada(controlcard,"SIGMAHEL",sigmahel,30.0d0)
2296 call reada(controlcard,"SIGMABET",sigmabet,40.0d0)
2297 call reada(controlcard,"WDIHC",wdihc,0.591d0)
2298 write (iout,*) "Weight of the dihedral restraint term",wdihc
2299 read(inp,'(9x,3f7.3)')&
2300 (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
2301 write (iout,*) "The secprob array"
2303 write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
2307 if (itype(i-3,molnum(i-3)).ne.ntyp1 .and. itype(i-2,molnum(i-2)).ne.ntyp1&
2308 .and. itype(i-1,molnum(i-1)).ne.ntyp1 .and. itype(i,molnum(i)).ne.ntyp1) then
2309 ndih_constr=ndih_constr+1
2310 idih_constr(ndih_constr)=i
2313 vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
2314 sumv=sumv+vpsipred(j,ndih_constr)
2317 vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
2319 phibound(1,ndih_constr)=phihel*deg2rad
2320 phibound(2,ndih_constr)=phibet*deg2rad
2321 sdihed(1,ndih_constr)=sigmahel*deg2rad
2322 sdihed(2,ndih_constr)=sigmabet*deg2rad
2326 'There are',ndih_constr,&
2327 ' bimodal restraints on gamma angles.'
2329 write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,&
2330 restyp(itype(idih_constr(i)-2,molnum(idih_constr(i)-2)),&
2331 molnum(idih_constr(i)-2)),idih_constr(i)-2,&
2332 restyp(itype(idih_constr(i)-1,molnum(idih_constr(i)-1)),&
2333 molnum(idih_constr(i)-2)),idih_constr(i)-1,&
2334 phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,&
2335 phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,&
2336 (vpsipred(j,i),j=1,3)
2339 endif ! endif ndif_constr.gt.0
2340 endif ! with_dihed_constr
2341 if (with_theta_constr) then
2342 !C with_theta_constr is keyword allowing for occurance of theta constrains
2343 read (inp,*) ntheta_constr
2344 !C ntheta_constr is the number of theta constrains
2345 if (ntheta_constr.gt.0) then
2346 !C read (inp,*) ftors
2347 read (inp,*) (itheta_constr(i),theta_constr0(i),&
2348 theta_drange(i),for_thet_constr(i),&
2350 !C the above code reads from 1 to ntheta_constr
2351 !C itheta_constr(i) residue i for which is theta_constr
2352 !C theta_constr0 the global minimum value
2353 !C theta_drange is range for which there is no energy penalty
2354 !C for_thet_constr is the force constant for quartic energy penalty
2357 'There are',ntheta_constr,' constraints on phi angles.'
2358 do i=1,ntheta_constr
2359 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),&
2364 do i=1,ntheta_constr
2365 theta_constr0(i)=deg2rad*theta_constr0(i)
2366 theta_drange(i)=deg2rad*theta_drange(i)
2368 !C if(me.eq.king.or..not.out1file)
2369 !C & write (iout,*) 'FTORS',ftors
2370 !C do i=1,ntheta_constr
2371 !C ii = itheta_constr(i)
2372 !C thetabound(1,ii) = phi0(i)-drange(i)
2373 !C thetabound(2,ii) = phi0(i)+drange(i)
2375 endif ! ntheta_constr.gt.0
2376 endif! with_theta_constr
2377 if (constr_homology.gt.0) then
2378 !c write (iout,*) "About to call read_constr_homology"
2380 call read_constr_homology
2381 !c write (iout,*) "Exit read_constr_homology"
2383 if (indpdb.gt.0 .or. pdbref) then
2387 c(j,i)=crefjlee(j,i)
2388 cref(j,i,kkk)=crefjlee(j,i)
2393 write (iout,*) "Array C"
2395 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),&
2398 write (iout,*) "Array Cref"
2400 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),&
2401 (cref(j,i+nres),j=1,3)
2405 call int_from_cart1(.false.)
2406 call sc_loc_geom(.false.)
2408 thetaref(i)=theta(i)
2410 write (iout,*) i," phiref",phiref(i)," thetaref",thetaref(i)
2414 dc(j,i)=c(j,i+1)-c(j,i)
2415 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
2420 dc(j,i+nres)=c(j,i+nres)-c(j,i)
2421 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
2430 write(iout,*)"przed ini_int_tab"
2432 write(iout,*)"po ini_int_tab"
2433 write(iout,*)"przed setup var"
2435 write(iout,*)"po setup var"
2436 if (ns.gt.0.and.dyn_ss) then
2440 forcon(i-nss)=forcon(i)
2447 ! call hpb_partition
2449 dyn_ss_mask(iss(i))=.true.
2453 write (iout,*) "molread: REFSTR",refstr
2455 if (.not.pdbref) then
2456 call read_angles(inp,*38)
2458 38 write (iout,'(a)') 'Error reading reference structure.'
2460 call mp_stopall(Error_Msg)
2462 stop 'Error reading reference structure'
2471 cref(j,i,kkk)=c(j,i)
2475 call contact(.true.,ncont_ref,icont_ref)
2478 end subroutine molread
2479 !-----------------------------------------------------------------------------
2480 subroutine openunits
2482 ! include 'DIMENSIONS'
2483 use control_data, only: from_cx,from_bx,from_cart
2487 character(len=3) :: liczba
2488 ! include "COMMON.MPI"
2490 ! include 'COMMON.IOUNITS'
2491 ! include 'COMMON.CONTROL'
2492 integer :: lenpre,lenpot !,ilen
2494 character(len=16) :: cformat,cprint
2495 ! character(len=16) ucase
2496 integer :: lenint,lenout
2497 call getenv('INPUT',prefix)
2498 call getenv('OUTPUT',prefout)
2499 call getenv('INTIN',prefintin)
2500 call getenv('COORD',cformat)
2501 call getenv('PRINTCOOR',cprint)
2502 call getenv('SCRATCHDIR',scratchdir)
2505 if (index(ucase(cformat),'CX').gt.0) then
2511 lenout=ilen(prefout)
2512 lenint=ilen(prefintin)
2513 ! Get the names and open the input files
2514 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
2516 write (liczba,'(bz,i3.3)') me
2517 outname=prefout(:lenout)//'_clust.out_'//liczba
2519 outname=prefout(:lenout)//'_clust.out'
2522 intinname=prefintin(:lenint)//'.bx'
2523 else if (from_cx) then
2524 intinname=prefintin(:lenint)//'.cx'
2526 intinname=prefintin(:lenint)//'.int'
2528 rmsname=prefintin(:lenint)//'.rms'
2529 open (jplot,file=prefout(:ilen(prefout))//'.tex',&
2531 open (jrms,file=rmsname,status='unknown')
2532 open(iout,file=outname,status='unknown')
2533 ! Get parameter filenames and open the parameter files.
2534 call getenv('BONDPAR',bondname)
2535 open (ibond,file=bondname,status='old')
2536 call getenv('THETPAR',thetname)
2537 open (ithep,file=thetname,status='old')
2538 call getenv('ROTPAR',rotname)
2539 open (irotam,file=rotname,status='old')
2540 call getenv('TORPAR',torname)
2541 open (itorp,file=torname,status='old')
2542 call getenv('TORDPAR',tordname)
2543 open (itordp,file=tordname,status='old')
2544 call getenv('FOURIER',fouriername)
2545 open (ifourier,file=fouriername,status='old')
2546 call getenv('ELEPAR',elename)
2547 open (ielep,file=elename,status='old')
2548 call getenv('SIDEPAR',sidename)
2549 open (isidep,file=sidename,status='old')
2550 call getenv('SIDEP',sidepname)
2551 open (isidep1,file=sidepname,status="old")
2552 call getenv('SCCORPAR',sccorname)
2553 open (isccor,file=sccorname,status="old")
2554 call getenv('BONDPAR_NUCL',bondname_nucl)
2555 open (ibond_nucl,file=bondname_nucl,status='old')
2556 call getenv('THETPAR_NUCL',thetname_nucl)
2557 open (ithep_nucl,file=thetname_nucl,status='old')
2558 call getenv('ROTPAR_NUCL',rotname_nucl)
2559 open (irotam_nucl,file=rotname_nucl,status='old')
2560 call getenv('TORPAR_NUCL',torname_nucl)
2561 open (itorp_nucl,file=torname_nucl,status='old')
2562 call getenv('TORDPAR_NUCL',tordname_nucl)
2563 open (itordp_nucl,file=tordname_nucl,status='old')
2564 call getenv('SIDEPAR_NUCL',sidename_nucl)
2565 open (isidep_nucl,file=sidename_nucl,status='old')
2566 call getenv('SIDEPAR_SCBASE',sidename_scbase)
2567 open (isidep_scbase,file=sidename_scbase,status='old')
2568 call getenv('PEPPAR_PEPBASE',pepname_pepbase)
2569 open (isidep_pepbase,file=pepname_pepbase,status='old')
2570 call getenv('SCPAR_PHOSPH',pepname_scpho)
2571 open (isidep_scpho,file=pepname_scpho,status='old')
2572 call getenv('PEPPAR_PHOSPH',pepname_peppho)
2573 open (isidep_peppho,file=pepname_peppho,status='old')
2576 call getenv('LIPTRANPAR',liptranname)
2577 open (iliptranpar,file=liptranname,status='old')
2578 call getenv('TUBEPAR',tubename)
2579 open (itube,file=tubename,status='old')
2580 call getenv('IONPAR',ionname)
2581 open (iion,file=ionname,status='old')
2582 call getenv('IONPAR_NUCL',ionnuclname)
2583 open (iionnucl,file=ionnuclname,status='old')
2584 call getenv('IONPAR_TRAN',iontranname)
2585 open (iiontran,file=iontranname,status='old')
2590 ! 8/9/01 In the newest version SCp interaction constants are read from a file
2591 ! Use -DOLDSCP to use hard-coded constants instead.
2593 call getenv('SCPPAR',scpname)
2594 open (iscpp,file=scpname,status='old')
2595 call getenv('SCPPAR_NUCL',scpname_nucl)
2596 open (iscpp_nucl,file=scpname_nucl,status='old')
2600 end subroutine openunits
2601 !-----------------------------------------------------------------------------
2603 !-----------------------------------------------------------------------------
2604 subroutine pdboutC(etot,rmsd,tytul)
2606 use energy_data, only: ihpb,jhpb,itype,molnum
2607 use energy, only:boxshift,to_box
2608 ! implicit real*8 (a-h,o-z)
2609 ! include 'DIMENSIONS'
2610 ! include 'COMMON.CONTROL'
2611 ! include 'COMMON.CHAIN'
2612 ! include 'COMMON.INTERACT'
2613 ! include 'COMMON.NAMES'
2614 ! include 'COMMON.IOUNITS'
2615 ! include 'COMMON.HEADER'
2616 ! include 'COMMON.SBRIDGE'
2617 ! include 'COMMON.TEMPFAC'
2618 character(len=50) :: tytul
2619 character(len=1) :: chainid(10)=(/'A','B','C','D','E','F',&
2621 integer :: ica(nres),k,iti1,ki
2622 real(kind=8) :: etot,rmsd
2623 integer :: iatom,ichain,ires,i,j,iti,mnum,mnum1,boxxxshift(3)
2624 real(kind=8) :: boxxxx(3),cbeg(3),Rdist(3)
2625 write (ipdb,'(3a,1pe15.5,a,0pf7.2)') 'REMARK ',tytul(:20),&
2626 ' ENERGY ',etot,' RMS ',rmsd
2630 ! boxxxshift(1)=int(c(1,nnt)/boxxsize)
2631 ! if (c(1,nnt).lt.0) boxxxshift(1)=boxxxshift(1)-1
2632 ! boxxxshift(2)=int(c(2,nnt)/boxzsize)
2633 ! if (c(2,nnt).lt.0) boxxxshift(2)=boxxxshift(2)-1
2634 ! boxxxshift(3)=int(c(3,nnt)/boxzsize)
2635 ! if (c(3,nnt).lt.0) boxxxshift(3)=boxxxshift(3)-1
2637 ! if (boxxxshift(i).gt.2) boxxxshift(i)=2
2638 ! if (boxxxshift(i).lt.-2) boxxxshift(i)=-2
2644 call to_box(cbeg(1),cbeg(2),cbeg(3))
2653 iti1=itype(i+1,mnum)
2655 if ((iti.eq.ntyp1_molec(mnum)).and.(iti1.eq.ntyp1_molec(mnum1))) cycle
2657 if (iti.eq.ntyp1_molec(mnum)) then
2660 write (ipdb,'(a)') 'TER'
2663 ! if (molnum(i).ge.1) then
2667 if (itype(i,mnum).ne.ntyp1_molec(mnum)) then
2673 if (itype(ki,mnum).eq.ntyp1_molec(mnum)) ki=ki+1
2674 if (itype(ki,mnum).eq.ntyp1_molec(mnum)) ki=ki+1
2680 call to_box(c(1,i),c(2,i),c(3,i))
2683 Rdist(k) = boxshift(c(k,i) - cbeg(k),boxxxx(k))
2684 c(k,i)=cbeg(k)+Rdist(k)
2686 if ((itype(i,molnum(i)).ne.10).and.(molnum(i).ne.5)) then
2687 call to_box(c(1,i+nres),c(2,i+nres),c(3,i+nres))
2689 Rdist(k) = boxshift(c(k,i+nres) - cbeg(k),boxxxx(k))
2690 c(k,i+nres)=cbeg(k)+Rdist(k)
2693 write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),&
2694 ires,(c(j,i),j=1,3),1.0d0,tempfac(1,i)
2695 if ((iti.ne.10).and.(mnum.ne.5)) then
2697 write (ipdb,20) iatom,restyp(iti,mnum),chainid(ichain),&
2698 ires,(c(j,nres+i),j=1,3),1.0d0,tempfac(2,i)
2702 write (ipdb,'(a)') 'TER'
2706 if ((itype(i,mnum).eq.ntyp1).or.(mnum.eq.5)) cycle
2707 if (itype(i,mnum).eq.10 .and. itype(i+1,mnum1).ne.ntyp1_molec(mnum1)) then
2708 write (ipdb,30) ica(i),ica(i+1)
2709 else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum1).ne.ntyp1_molec(mnum1)) then
2710 write (ipdb,30) ica(i),ica(i+1),ica(i)+1
2711 else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) then
2712 write (ipdb,30) ica(i),ica(i)+1
2715 if ((itype(nct,molnum(nct)).ne.10).and.(molnum(nct).ne.5)) then
2716 write (ipdb,30) ica(nct),ica(nct)+1
2719 write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
2721 write (ipdb,'(a6)') 'ENDMDL'
2722 10 FORMAT ('ATOM',I7,' CA ',A3,1X,A1,I4,4X,3F8.3,2f6.2)
2723 20 FORMAT ('ATOM',I7,' CB ',A3,1X,A1,I4,4X,3F8.3,2f6.2)
2724 30 FORMAT ('CONECT',8I5)
2726 end subroutine pdboutC
2727 !-----------------------------------------------------------------------------
2728 subroutine cartout(igr,i,etot,free,rmsd,plik)
2729 ! implicit real*8 (a-h,o-z)
2730 ! include 'DIMENSIONS'
2731 ! include 'sizesclu.dat'
2732 ! include 'COMMON.IOUNITS'
2733 ! include 'COMMON.CHAIN'
2734 ! include 'COMMON.VAR'
2735 ! include 'COMMON.LOCAL'
2736 ! include 'COMMON.INTERACT'
2737 ! include 'COMMON.NAMES'
2738 ! include 'COMMON.GEO'
2739 ! include 'COMMON.CLUSTER'
2740 integer :: igr,i,j,k
2741 real(kind=8) :: etot,free,rmsd
2742 character(len=80) :: plik
2743 open (igeom,file=plik,position='append')
2744 write (igeom,'(2e15.5,f10.5,$)') etot,free,rmsd
2745 write (igeom,'(i4,$)') &
2746 nss_all(i),(ihpb_all(j,i),jhpb_all(j,i),j=1,nss_all(i))
2747 write (igeom,'(i10)') iscore(i)
2748 write (igeom,'(8f10.5)') &
2749 ((allcart(k,j,i),k=1,3),j=1,nres),&
2750 ((allcart(k,j+nres,i),k=1,3),j=nnt,nct)
2752 end subroutine cartout
2753 !------------------------------------------------------------------------------
2754 subroutine read_klapaucjusz
2756 use control_data, only:unres_pdb
2757 use geometry_data, only:theta,thetaref,phi,phiref,&
2760 ! include 'DIMENSIONS'
2764 ! include 'COMMON.SETUP'
2765 ! include 'COMMON.CONTROL'
2766 ! include 'COMMON.HOMOLOGY'
2767 ! include 'COMMON.CHAIN'
2768 ! include 'COMMON.IOUNITS'
2769 ! include 'COMMON.MD'
2770 ! include 'COMMON.GEO'
2771 ! include 'COMMON.INTERACT'
2772 ! include 'COMMON.NAMES'
2773 character(len=256) fragfile
2774 integer, dimension(:), allocatable :: ninclust,nresclust,itype_temp,&
2776 integer, dimension(:,:), allocatable :: iresclust,inclust
2779 character(len=2) :: kic2
2780 character(len=24) :: model_ki_dist, model_ki_angle
2781 character(len=500) :: controlcard
2782 integer :: ki, i, j, jj,k, l, i_tmp,&
2784 ik,ll,lll,ii_old,ii,iii,ichain,kk,iistart,iishift,lim_xx,igr,&
2785 i01,i10,nnt_chain,nct_chain
2786 real(kind=8) :: distal
2787 logical :: lprn = .true.
2788 integer :: nres_temp
2791 logical :: liiflag,lfirst
2793 real(kind=8) rescore_tmp,x12,y12,z12,rescore2_tmp,dist_cut
2794 real(kind=8), dimension (:,:), allocatable :: rescore
2795 real(kind=8), dimension (:,:), allocatable :: rescore2
2796 character(len=24) :: tpl_k_rescore
2797 character(len=256) :: pdbfile
2800 ! For new homol impl
2802 ! include 'COMMON.VAR'
2804 ! write (iout,*) "READ_KLAPAUCJUSZ"
2805 ! print *,"READ_KLAPAUCJUSZ"
2807 call getenv("FRAGFILE",fragfile)
2808 write (iout,*) "Opening", fragfile
2810 open(ientin,file=fragfile,status="old",err=10)
2811 ! write (iout,*) " opened"
2820 itype_temp(:)=itype(:,1)
2824 ! write (iout,*) "Entering loop"
2827 DO IGR = 1,NCHAIN_GROUP
2829 ! write (iout,*) "igr",igr
2831 read(ientin,*) constr_homology,nclust
2832 if (start_from_model) then
2833 nmodel_start=constr_homology
2840 ichain=iequiv(1,igr)
2841 nnt_chain=chain_border(1,ichain)-chain_border1(1,ichain)+1
2842 nct_chain=chain_border(2,ichain)-chain_border1(1,ichain)+1
2843 ! write (iout,*) "nnt_chain",nnt_chain," nct_chain",nct_chain
2846 if(.not.allocated(pdbfiles_chomo)) allocate(pdbfiles_chomo(constr_homology))
2847 if(.not.allocated(nres_chomo)) allocate(nres_chomo(constr_homology))
2848 do k=1,constr_homology
2849 read(ientin,'(a)') pdbfile
2850 write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file', &
2851 pdbfile(:ilen(pdbfile))
2852 pdbfiles_chomo(k)=pdbfile
2853 open(ipdbin,file=pdbfile,status='old',err=33)
2855 33 write (iout,'(a,5x,a)') 'Error opening PDB file',&
2856 pdbfile(:ilen(pdbfile))
2861 call readpdb_template(k)
2871 read(ientin,*) ninclust(i),nresclust(i)
2872 read(ientin,*) (inclust(k,i),k=1,ninclust(i))
2873 read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
2876 ! Loop over clusters
2879 do ll = 1,ninclust(l)
2882 ! write (iout,*) "l",l," ll",ll," k",k
2888 idomain(k,iresclust(i,l)+1) = 1
2890 idomain(k,iresclust(i,l)) = 1
2894 ! Distance restraints
2897 ! Copy the coordinates from reference coordinates (?)
2903 ! write (iout,*) "c(",j,i,") =",c(j,i)
2906 call int_from_cart(.true.,.false.)
2907 call sc_loc_geom(.false.)
2909 thetaref(i)=theta(i)
2913 if (waga_dist.ne.0.0d0) then
2916 do i = nnt_chain,nct_chain-2
2923 distal=dsqrt(x12*x12+y12*y12+z12*z12)
2924 ! write (iout,*) k,i,j,distal,dist2_cut
2926 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 &
2927 .and. distal.le.dist2_cut ) then
2933 ! write (iout,*) "k",k
2934 ! write (iout,*) "i",i," j",j," constr_homology",
2936 ires_homo(ii)=i+chain_border1(1,igr)-1
2937 jres_homo(ii)=j+chain_border1(1,igr)-1
2939 if (read2sigma) then
2942 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
2944 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
2945 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = &
2946 sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2948 if (odl(k,ii).le.dist_cut) then
2949 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
2952 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
2953 dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
2955 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
2956 dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2960 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
2963 ! l_homo(k,ii)=.false.
2970 ! Theta, dihedral and SC retraints
2972 if (waga_angle.gt.0.0d0) then
2973 do i = nnt_chain+3,nct_chain
2974 iii=i+chain_border1(1,igr)-1
2975 if (idomain(k,i).eq.0) then
2976 ! sigma_dih(k,i)=0.0
2979 dih(k,iii)=phiref(i)
2981 (rescore(k,i)+rescore(k,i-1)+ &
2982 rescore(k,i-2)+rescore(k,i-3))/4.0
2983 ! write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
2984 ! & " sigma_dihed",sigma_dih(k,i)
2985 if (sigma_dih(k,iii).ne.0) &
2986 sigma_dih(k,iii)=1.0d0/(sigma_dih(k,iii)*sigma_dih(k,iii))
2991 if (waga_theta.gt.0.0d0) then
2992 do i = nnt_chain+2,nct_chain
2993 iii=i+chain_border1(1,igr)-1
2994 if (idomain(k,i).eq.0) then
2995 ! sigma_theta(k,i)=0.0
2998 thetatpl(k,iii)=thetaref(i)
2999 sigma_theta(k,iii)=(rescore(k,i)+rescore(k,i-1)+ &
3001 if (sigma_theta(k,iii).ne.0) &
3002 sigma_theta(k,iii)=1.0d0/ &
3003 (sigma_theta(k,iii)*sigma_theta(k,iii))
3007 if (waga_d.gt.0.0d0) then
3008 do i = nnt_chain,nct_chain
3009 iii=i+chain_border1(1,igr)-1
3010 if (itype(i,1).eq.10) cycle
3011 if (idomain(k,i).eq.0 ) then
3015 xxtpl(k,iii)=xxref(i)
3016 yytpl(k,iii)=yyref(i)
3017 zztpl(k,iii)=zzref(i)
3018 sigma_d(k,iii)=rescore(k,i)
3019 if (sigma_d(k,iii).ne.0) &
3020 sigma_d(k,iii)=1.0d0/(sigma_d(k,iii)*sigma_d(k,iii))
3021 ! if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
3027 ! remove distance restraints not used in any model from the list
3028 ! shift data in all arrays
3030 ! write (iout,*) "ii_old",ii_old
3031 if (waga_dist.ne.0.0d0) then
3033 write (iout,*) "Distance restraints from templates"
3035 write(iout,'(4i5,100(2f8.2,1x,l1,4x))') &
3036 iii,ii_in_use(iii),ires_homo(iii),jres_homo(iii), &
3037 (odl(ki,iii),1.0d0/dsqrt(sigma_odl(ki,iii)),l_homo(ki,iii), &
3038 ki=1,constr_homology)
3044 do i=nnt_chain,nct_chain-2
3047 ! if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3048 ! & .and. distal.le.dist2_cut ) then
3049 ! write (iout,*) "i",i," j",j," ii",ii," i_in_use",ii_in_use(ii)
3051 if (ii_in_use(ii).eq.0.and.liiflag.or. &
3052 ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
3059 if(i10.eq.lim_odl) i10=i10+1
3061 ires_homo(iistart+ki)=ires_homo(ki+i01)
3062 jres_homo(iistart+ki)=jres_homo(ki+i01)
3063 ii_in_use(iistart+ki)=ii_in_use(ki+i01)
3064 do k=1,constr_homology
3065 odl(k,iistart+ki)=odl(k,ki+i01)
3066 sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
3067 l_homo(k,iistart+ki)=l_homo(k,ki+i01)
3070 iistart=iistart+i10-i01
3073 if (ii_in_use(ii).ne.0.and..not.liiflag) then
3086 ichain=iequiv(i,igr)
3088 do j=nnt_chain,nct_chain
3089 jj=j+chain_border1(1,ichain)-chain_border1(1,iequiv(1,igr))
3090 do k=1,constr_homology
3092 sigma_dih(k,jj)=sigma_dih(k,j)
3093 thetatpl(k,jj)=thetatpl(k,j)
3094 sigma_theta(k,jj)=sigma_theta(k,j)
3095 xxtpl(k,jj)=xxtpl(k,j)
3096 yytpl(k,jj)=yytpl(k,j)
3097 zztpl(k,jj)=zztpl(k,j)
3098 sigma_d(k,jj)=sigma_d(k,j)
3102 jj=chain_border1(1,ichain)-chain_border1(1,iequiv(i-1,igr))
3103 ! write (iout,*) "igr",igr," i",i," ichain",ichain," jj",jj
3104 do j=ii_old+1,lim_odl
3105 ires_homo(j+lll)=ires_homo(j)+jj
3106 jres_homo(j+lll)=jres_homo(j)+jj
3107 do k=1,constr_homology
3108 odl(k,j+lll)=odl(k,j)
3109 sigma_odl(k,j+lll)=sigma_odl(k,j)
3110 l_homo(k,j+lll)=l_homo(k,j)
3121 if (waga_angle.gt.0.0d0) lim_dih=nct-nnt-2
3123 itype(:,1)=itype_temp(:)
3126 10 stop "Error in fragment file"
3127 end subroutine read_klapaucjusz
3129 !-----------------------------------------------------------------------------