2 !-----------------------------------------------------------------------------
6 use io_base !, only: ilen
7 use geometry_data, only: nres,c
8 use energy_data, only: nnt,nct,nss
9 use control_data, only: lside
11 !-----------------------------------------------------------------------------
14 !-----------------------------------------------------------------------------
16 !-----------------------------------------------------------------------------
18 !-----------------------------------------------------------------------------
19 SUBROUTINE WRTCLUST(NCON,ICUT,PRINTANG,PRINTPDB,printmol2,ib)
21 use hc_, only: ioffset
22 use control_data, only: lprint_cart,lprint_int,titel
23 use geometry, only: int_from_cart1,nres
24 ! implicit real*8 (a-h,o-z)
25 ! include 'DIMENSIONS'
26 ! include 'sizesclu.dat'
27 integer,parameter :: num_in_line=5
28 LOGICAL :: PRINTANG(max_cut)
29 integer :: PRINTPDB(max_cut),printmol2(max_cut)
30 ! include 'COMMON.CONTROL'
31 ! include 'COMMON.HEADER'
32 ! include 'COMMON.CHAIN'
33 ! include 'COMMON.VAR'
34 ! include 'COMMON.CLUSTER'
35 ! include 'COMMON.IOUNITS'
36 ! include 'COMMON.GEO'
37 ! include 'COMMON.FREE'
38 ! include 'COMMON.TEMPFAC'
39 real(kind=8) :: rmsave(maxgr)
40 CHARACTER(len=64) :: prefixp,NUMM,MUMM,EXTEN,extmol
41 character(len=80) :: cfname
42 character(len=8) :: ctemper
43 DATA EXTEN /'.pdb'/,extmol /'.mol2'/,NUMM /'000'/,&
46 integer :: ncon,icut,ib
47 integer :: i,ii,ii1,ii2,igr,ind1,ind2,ico,icon,&
48 irecord,nrecord,j,k,jj,ind,ncon_lim,ncon_out
49 real(kind=8) :: temper,curr_dist,emin,qpart,boltz,&
50 ave_dim,amax_dim,emin1
53 allocate(tempfac(2,nres))
58 ! print *,"calling WRTCLUST",ncon
59 ! write (iout,*) "ICUT",icut," PRINTPDB ",PRINTPDB(icut)
62 temper=1.0d0/(beta_h(ib)*1.987d-3)
63 if (temper.lt.100.0d0) then
64 write(ctemper,'(f3.0)') temper
66 else if (temper.lt.1000.0) then
67 write (ctemper,'(f4.0)') temper
70 write (ctemper,'(f5.0)') temper
74 do i=1,ncon*(ncon-1)/2
77 close(80,status='delete')
79 ! PRINT OUT THE RESULTS OF CLUSTER ANALYSIS
81 ii1= index(intinname,'/')
86 ii2=index(intinname(ii1:),'/')
88 ii = ii1+index(intinname(ii1:),'.')-1
94 prefixp=intinname(ii1:ii)
95 !d print *,icut,printang(icut),printpdb(icut),printmol2(icut)
96 !d print *,'ecut=',ecut
99 WRITE (iout,200) IGR,totfree_gr(igr)/beta_h(ib),LICZ(IGR)
100 NRECORD=LICZ(IGR)/num_in_line
102 DO 63 IRECORD=1,NRECORD
103 IND2=IND1+num_in_line-1
104 WRITE (iout,300) (list_conf(NCONF(IGR,ICO)),&
105 totfree(NCONF(IGR,ICO))/beta_h(ib),ICO=IND1,IND2)
108 WRITE (iout,300) (list_conf(NCONF(IGR,ICO)),&
109 totfree(NCONF(IGR,ICO))/beta_h(ib),ICO=IND1,LICZ(IGR))
111 ICON=list_conf(NCONF(IGR,1))
112 ! WRITE (iout,'(16F5.0)') (rad2deg*phiall(IND,icon),IND=4,nphi+3)
113 ! 12/8/93 Estimation of "diameters" of the subsequent families.
116 ! write (iout,*) "ecut",ecut
119 if (totfree(ii)-emin .gt. ecut) goto 10
124 ind=ioffset(ncon,ii,jj)
126 ind=ioffset(ncon,jj,ii)
128 ! write (iout,*) " ncon",ncon,"i",i," j",j," ii",ii," jj",jj,
131 curr_dist=dabs(diss(ind)+0.0d0)
132 ! write(iout,'(i10,4i4,f12.4)') ind,ii,jj,list_conf(ii),
133 ! & list_conf(jj),curr_dist
134 if (curr_dist .gt. amax_dim) amax_dim=curr_dist
135 ave_dim=ave_dim+curr_dist**2
138 10 if (licz(igr) .gt. 1) &
139 ave_dim=sqrt(ave_dim/(licz(igr)*(licz(igr)-1)/2))
140 write (iout,'(/A,F8.1,A,F8.1)') &
141 'Max. distance in the family:',amax_dim,&
142 '; average distance in the family:',ave_dim
147 boltz=dexp(-totfree(icon))
148 rmsave(igr)=rmsave(igr)+boltz*rmstb(icon)
151 rmsave(igr)=rmsave(igr)/qpart
152 write (iout,'(a,f5.2,a)') "Average RMSD",rmsave(igr)," A"
155 WRITE (iout,500) (list_conf(I),IASS(I),I=1,NCON)
156 ! print *,icut,printang(icut)
157 IF (PRINTANG(ICUT) .and. (lprint_cart .or. lprint_int)) then
159 ! print *,'emin',emin,' ngr',ngr
160 if (lprint_cart) then
161 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))&
164 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))&
169 if (totfree_gr(igr)-emin.le.ecut) then
170 if (lprint_cart) then
171 call cartout(igr,icon,totfree(icon)/beta_h(ib),&
172 totfree_gr(igr)/beta_h(ib),&
175 ! print '(a)','calling briefout'
178 c(j,i)=allcart(j,i,icon)
181 call int_from_cart1(.false.)
182 !el call briefout(igr,iscore(icon),totfree(icon)/beta_h(ib),&
183 !el totfree_gr(igr),nss_all(icon),ihpb_all(1,icon),&
184 !el jhpb_all(1,icon),cfname)
185 call briefout(igr,totfree(icon)/beta_h(ib),&
187 ! print '(a)','exit briefout'
193 IF (PRINTPDB(ICUT).gt.0) THEN
194 ! Write out a number of conformations from each family in PDB format and
195 ! create InsightII command file for their displaying in different colors
196 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))&
198 write (iout,*) "cfname",cfname
199 OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
200 write (ipdb,'(a,f8.2)') &
201 "REMAR AVERAGE CONFORMATIONS AT TEMPERATURE",temper
207 DO WHILE(I.LE.NGR .AND. totfree_gr(i)-EMIN.LE.ECUT)
208 ! write (iout,*) "i",i," ngr",ngr,totfree_gr(I),EMIN,ecut
209 write (NUMM,'(bz,i4.4)') i
210 ncon_lim=min0(licz(i),printpdb(icut))
211 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))&
212 //"K_"//numm(:ilen(numm))//exten
213 OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
214 write (ipdb,'("REMARK CLUSTER",i5," FREE ENERGY",1pe14.5," AVE RMSD",0pf5.2)')&
215 i,totfree_gr(i)/beta_h(ib),rmsave(i) !'
216 ! Write conformations of the family i to PDB files
218 do while (ncon_out.lt.printpdb(icut) .and. &
219 ncon_out.lt.licz(i).and. &
220 totfree(nconf(i,ncon_out+1))-EMIN1.LE.ECUT)
222 ! write (iout,*) i,ncon_out,nconf(i,ncon_out),
223 ! & totfree(nconf(i,ncon_out)),emin1,ecut
225 write (iout,*) "ncon_out",ncon_out
235 c(k,ii)=allcart(k,ii,icon)
238 call pdboutC(totfree(icon)/beta_h(ib),rmstb(icon),titel)
239 write (ipdb,'("TER")')
242 ! Average structures and structures closest to average
243 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))&
245 OPEN(ipdb,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED',&
248 write (ipdb,'(a,i5)') "REMARK CLUSTER",i
249 call pdboutC(totfree_gr(i)/beta_h(ib),rmsave(i),titel)
250 write (ipdb,'("TER")')
251 call closest_coord(i)
252 call pdboutC(totfree_gr(i)/beta_h(ib),rmsave(i),titel)
253 write (ipdb,'("TER")')
260 IF (printmol2(icut).gt.0) THEN
261 ! Write out a number of conformations from each family in PDB format and
262 ! create InsightII command file for their displaying in different colors
267 DO WHILE(I.LE.NGR .AND. totfree_gr(i)-EMIN.LE.ECUT)
268 write (NUMM,'(bz,i4.4)') i
269 cfname=prefixp(:ilen(prefixp))//"_T"//ctemper(:ilen(ctemper))&
270 //"K_"//numm(:ilen(numm))//extmol
271 OPEN(imol2,FILE=CFNAME,STATUS='UNKNOWN',FORM='FORMATTED')
273 do while (ncon_out.lt.printmol2(icut) .and. &
274 ncon_out.lt.licz(i).and. &
275 totfree(nconf(i,ncon_out+1))-EMIN1.LE.ECUT)
282 c(k,ii)=allcart(k,ii,icon)
285 CALL MOL2OUT(totfree(icon)/beta_h(ib),'STRUCTURE'//numm)
293 100 FORMAT (//'THERE ARE ',I4,' FAMILIES OF CONFORMATIONS')
294 200 FORMAT (/'FAMILY ',I4,' WITH TOTAL FREE ENERGY',1pE15.5,&
295 ' CONTAINS ',I4,' CONFORMATION(S): ')
296 ! 300 FORMAT ( 8(I4,F6.1))
297 300 FORMAT (5(I4,1pe12.3))
298 400 FORMAT (//'ASSIGNMENT OF CONSECUTIVE CONFORMATIONS TO FAMILIES:')
299 500 FORMAT (8(2I4,2X))
300 600 FORMAT ('REMARK FAMILY',I4,' CONFORMATION',I4,' ENERGY ',E15.6)
302 END SUBROUTINE WRTCLUST
303 !------------------------------------------------------------------------------
304 subroutine ave_coord(igr)
306 use control_data, only:lside
307 use regularize_, only:fitsq,matvec
309 ! include 'DIMENSIONS'
310 ! include 'sizesclu.dat'
311 ! include 'COMMON.CONTROL'
312 ! include 'COMMON.CLUSTER'
313 ! include 'COMMON.CHAIN'
314 ! include 'COMMON.INTERACT'
315 ! include 'COMMON.VAR'
316 ! include 'COMMON.TEMPFAC'
317 ! include 'COMMON.IOUNITS'
319 real(kind=8) :: przes(3),obrot(3,3)
320 real(kind=8) :: xx(3,2*nres),yy(3,2*nres),csq(3,2*nres) !(3,maxres2)
322 integer :: i,ii,j,k,icon,jcon,igr
323 real(kind=8) :: rms,boltz,qpart,cwork(3,2*nres),cref1(3,2*nres) !(3,maxres2)
324 ! write (iout,*) "AVE_COORD: igr",igr
327 boltz = dexp(-totfree(jcon)+eref)
331 c(j,i)=allcart(j,i,jcon)*boltz
332 cref1(j,i)=allcart(j,i,jcon)
333 csq(j,i)=allcart(j,i,jcon)**2*boltz
343 xx(j,ii)=allcart(j,i,jcon)
348 ! if (itype(i).ne.10) then
351 xx(j,ii)=allcart(j,i+nres,jcon)
352 yy(j,ii)=cref1(j,i+nres)
356 call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
360 cwork(j,i)=allcart(j,i,jcon)
363 call fitsq(rms,cwork(1,nnt),cref1(1,nnt),nct-nnt+1,przes,obrot &
366 ! write (iout,*) "rms",rms
368 ! write (iout,'(i3,f10.5,5x,3f10.5)')i,przes(i),(obrot(i,j),j=1,3)
371 print *,'error, rms^2 = ',rms,icon,jcon
374 if (non_conv) print *,non_conv,icon,jcon
375 boltz=dexp(-totfree(jcon)+eref)
376 qpart = qpart + boltz
379 xx(j,i)=allcart(j,i,jcon)
382 call matvec(cwork,obrot,xx,2*nres)
384 ! write (iout,'(i5,2(3f10.5,5x))') i,(cwork(j,i),j=1,3),
385 ! & (allcart(j,i,jcon),j=1,3)
387 cwork(j,i)=cwork(j,i)+przes(j)
388 c(j,i)=c(j,i)+cwork(j,i)*boltz
389 csq(j,i)=csq(j,i)+cwork(j,i)**2*boltz
396 csq(j,i)=csq(j,i)/qpart-c(j,i)**2
398 ! write (iout,'(i5,3f10.5)') i,(csq(j,i),j=1,3)
404 tempfac(1,i)=tempfac(1,i)+csq(j,i)
405 tempfac(2,i)=tempfac(2,i)+csq(j,i+nres)
407 tempfac(1,i)=dsqrt(tempfac(1,i))
408 tempfac(2,i)=dsqrt(tempfac(2,i))
411 end subroutine ave_coord
412 !------------------------------------------------------------------------------
413 subroutine closest_coord(igr)
415 use regularize_, only:fitsq
417 ! include 'DIMENSIONS'
418 ! include 'sizesclu.dat'
419 ! include 'COMMON.IOUNITS'
420 ! include 'COMMON.CONTROL'
421 ! include 'COMMON.CLUSTER'
422 ! include 'COMMON.CHAIN'
423 ! include 'COMMON.INTERACT'
424 ! include 'COMMON.VAR'
426 real(kind=8) :: przes(3),obrot(3,3)
427 real(kind=8) :: xx(3,2*nres),yy(3,2*nres) !(3,maxres2)
428 integer :: i,ii,j,k,icon,jcon,jconmin,igr
429 real(kind=8) :: rms,rmsmin,cwork(3,2*nres)
439 xx(j,ii)=allcart(j,i,jcon)
444 ! if (itype(i).ne.10) then
447 xx(j,ii)=allcart(j,i+nres,jcon)
452 call fitsq(rms,xx(1,1),yy(1,1),ii,przes,obrot,non_conv)
456 cwork(j,i)=allcart(j,i,jcon)
459 call fitsq(rms,cwork(1,nnt),c(1,nnt),nct-nnt+1,przes,obrot,&
463 print *,'error, rms^2 = ',rms,icon,jcon
466 ! write (iout,*) "jcon",jcon," rms",rms," rmsmin",rmsmin
467 if (non_conv) print *,non_conv,icon,jcon
468 if (rms.lt.rmsmin) then
473 ! write (iout,*) "rmsmin",rmsmin," rms",rms
477 c(j,i)=allcart(j,i,jconmin)
481 end subroutine closest_coord
482 !-----------------------------------------------------------------------------
484 !-----------------------------------------------------------------------------
485 subroutine read_coords(ncon,*)
487 use energy_data, only: ihpb,jhpb,max_ene
488 use control_data, only: from_bx,from_cx
489 use control, only: tcpu
491 ! include "DIMENSIONS"
492 ! include "sizesclu.dat"
496 integer :: IERROR,ERRCODE !,STATUS(MPI_STATUS_SIZE)
497 ! include "COMMON.MPI"
499 use MPI_data, only: nprocs
501 ! include "COMMON.CONTROL"
502 ! include "COMMON.CHAIN"
503 ! include "COMMON.INTERACT"
504 ! include "COMMON.IOUNITS"
505 ! include "COMMON.VAR"
506 ! include "COMMON.SBRIDGE"
507 ! include "COMMON.GEO"
508 ! include "COMMON.CLUSTER"
509 character(len=3) :: liczba
511 integer :: i,j,jj,jjj,jj_old,icount,k,kk,l,ii,if,ib,&
513 integer :: ixdrf,iret,itmp
514 real(kind=4) :: prec,reini,refree,rmsdev
515 integer :: nrec,nlines,iscor,lenrec,lenrec_in
516 real(kind=8) :: energ,t_acq !,tcpu
517 !el integer ilen,iroof
518 !el external ilen,iroof
519 real(kind=8) :: rjunk
520 integer :: ntot_all(0:nprocs-1) !(0:maxprocs-1)
522 real(kind=8) :: energia(0:max_ene),etot
523 real(kind=4) :: csingle(3,2*nres+2)
524 integer :: Previous,Next
525 character(len=256) :: bprotfiles
526 ! print *,"Processor",me," calls read_protein_data"
528 if (me.eq.master) then
529 Previous=MPI_PROC_NULL
533 if (me.eq.nprocs-1) then
538 ! Set the scratchfile names
539 write (liczba,'(bz,i3.3)') me
541 allocate(STATUS(MPI_STATUS_SIZE))
543 ! 1/27/05 AL Change stored coordinates to single precision and don't store
544 ! energy components in the binary databases.
545 lenrec=12*(nres+nct-nnt+1)+4*(2*nss+2)+16
546 lenrec_in=12*(nres+nct-nnt+1)+4*(2*nss+2)+24
548 write (iout,*) "nres",nres," nnt",nnt," nct",nct," nss", nss
549 write (iout,*) "lenrec_in",lenrec_in
551 bprotfiles=scratchdir(:ilen(scratchdir))// &
552 "/"//prefix(:ilen(prefix))//liczba//".xbin"
554 ! allocate cluster arrays
555 allocate(energy(0:maxconf),totfree(0:maxconf)) !(0:maxconf)
556 allocate(entfac(maxconf)) !(maxconf)
557 allocate(rmstb(maxconf)) !(maxconf)
558 allocate(allcart(3,2*nres,maxstr_proc)) !(3,maxres2,maxstr_proc)
559 allocate(nss_all(maxstr_proc)) !(maxstr_proc)
560 allocate(ihpb_all(maxss,maxstr_proc),jhpb_all(maxss,maxstr_proc))!(maxss,maxstr_proc)
561 allocate(iscore(maxconf)) !(maxconf)
567 if (from_cart .and. .not. from_bx .and. .not. from_cx) then
569 read (intin,*,end=13,err=11) energy(icon),totfree(icon),&
571 nss_all(icon),(ihpb_all(ii,icon),jhpb_all(i,icon),&
572 i=1,nss_all(icon)),iscore(icon)
574 read (intin,*,end=13,err=11) energy(icon),rmstb(icon),&
575 nss_all(icon),(ihpb_all(ii,icon),jhpb_all(i,icon),&
576 i=1,nss_all(icon)),iscore(icon)
578 read (intin,'(8f10.5)',end=13,err=10) &
579 ((allcart(j,i,icon),j=1,3),i=1,nres),&
580 ((allcart(j,i+nres,icon),j=1,3),i=nnt,nct)
581 print *,icon,energy(icon),nss_all(icon),rmstb(icon)
583 read(intin,'(a80)',end=13,err=12) lineh
584 read(lineh(:5),*,err=8) ic
586 read(lineh(6:),*,err=8) energy(icon)
588 read(lineh(6:),*,err=8) energy(icon)
592 print *,'error, assuming e=1d10',lineh
596 !old read(lineh(18:),*,end=13,err=11) nss_all(icon)
597 ii = index(lineh(15:)," ")+15
598 read(lineh(ii:),*,end=13,err=11) nss_all(icon)
599 IF (NSS_all(icon).LT.9) THEN
600 read (lineh(20:),*,end=102) &
601 (IHPB_all(I,icon),JHPB_all(I,icon),I=1,NSS_all(icon)),&
604 read (lineh(20:),*,end=102) &
605 (IHPB_all(I,icon),JHPB_all(I,icon),I=1,8)
606 read (intin,*) (IHPB_all(I,icon),JHPB_all(I,icon),&
607 I=9,NSS_all(icon)),iscore(icon)
612 PRINT *,'IC:',IC,' ENERGY:',ENERGY(ICON)
613 call read_angles(intin,*13)
615 phiall(i,icon)=phi(i)
616 thetall(i,icon)=theta(i)
617 alphall(i,icon)=alph(i)
618 omall(i,icon)=omeg(i)
624 ! CALCULATE DISTANCES
626 10 print *,'something wrong with angles'
628 11 print *,'something wrong with NSS',nss
630 12 print *,'something wrong with header'
637 open (icbase,file=bprotfiles,status="unknown",&
638 form="unformatted",access="direct",recl=lenrec)
639 ! Read conformations from binary DA files (one per batch) and write them to
640 ! a binary DA scratchfile.
644 write (liczba,'(bz,i3.3)') me
645 IF (ME.EQ.MASTER) THEN
646 ! Only the master reads the database; it'll send it to the other procs
654 open (intin,file=intinname,status="old",form="unformatted",&
655 access="direct",recl=lenrec_in)
657 else if (from_cx) then
658 #if (defined(AIX) && !defined(JUBL))
659 call xdrfopen_(ixdrf,intinname, "r", iret)
661 call xdrfopen(ixdrf,intinname, "r", iret)
664 write (iout,*) "xdrfopen: iret",iret
666 write (iout,*) "Error: coordinate file ",&
667 intinname(:ilen(intinname))," does not exist."
670 call MPI_ABORT(MPI_COMM_WORLD,IERROR,ERRCODE)
675 write (iout,*) "Error: coordinate format not specified"
678 call MPI_ABORT(MPI_COMM_WORLD,IERROR,ERRCODE)
686 write (iout,*) "Opening file ",intinname(:ilen(intinname))
687 write (iout,*) "lenrec",lenrec_in
691 ! write (iout,*) "maxconf",maxconf
695 !el if (i.gt.maxconf) then
696 !el write (iout,*) "Error: too many conformations ",&
697 !el "(",maxconf,") maximum."
699 !el call MPI_Abort(MPI_COMM_WORLD,errcode,ierror)
703 ! write (iout,*) "i",i
706 read(intin,err=101,end=101) &
707 ((csingle(l,k),l=1,3),k=1,nres),&
708 ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
709 nss,(ihpb(k),jhpb(k),k=1,nss),&
711 entfac(jj+1),rmstb(jj+1),iscor
718 #if (defined(AIX) && !defined(JUBL))
719 call xdrf3dfcoord_(ixdrf, csingle, itmp, prec, iret)
720 if (iret.eq.0) goto 101
721 call xdrfint_(ixdrf, nss, iret)
722 if (iret.eq.0) goto 101
724 call xdrfint_(ixdrf, ihpb(j), iret)
725 if (iret.eq.0) goto 101
726 call xdrfint_(ixdrf, jhpb(j), iret)
727 if (iret.eq.0) goto 101
729 call xdrffloat_(ixdrf,reini,iret)
730 if (iret.eq.0) goto 101
731 call xdrffloat_(ixdrf,refree,iret)
732 if (iret.eq.0) goto 101
733 call xdrffloat_(ixdrf,rmsdev,iret)
734 if (iret.eq.0) goto 101
735 call xdrfint_(ixdrf,iscor,iret)
736 if (iret.eq.0) goto 101
738 ! write (iout,*) "calling xdrf3dfcoord"
739 call xdrf3dfcoord(ixdrf, csingle, itmp, prec, iret)
740 ! write (iout,*) "iret",iret
742 if (iret.eq.0) goto 101
743 call xdrfint(ixdrf, nss, iret)
744 ! write (iout,*) "iret",iret
745 ! write (iout,*) "nss",nss
747 if (iret.eq.0) goto 101
749 call xdrfint(ixdrf, ihpb(k), iret)
750 if (iret.eq.0) goto 101
751 call xdrfint(ixdrf, jhpb(k), iret)
752 if (iret.eq.0) goto 101
754 call xdrffloat(ixdrf,reini,iret)
755 if (iret.eq.0) goto 101
756 call xdrffloat(ixdrf,refree,iret)
757 if (iret.eq.0) goto 101
758 call xdrffloat(ixdrf,rmsdev,iret)
759 if (iret.eq.0) goto 101
760 call xdrfint(ixdrf,iscor,iret)
761 if (iret.eq.0) goto 101
773 c(l,nres+k)=csingle(l,nres+k-nnt+1)
778 write (iout,'(5hREAD ,i5,3f15.4,i10)') &
779 jj+1,energy(jj+1),entfac(jj+1),&
781 write (iout,*) "Conformation",jjj+1,jj+1
782 write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
783 write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
786 call add_new_cconf(jjj,jj,jj_old,icount,Next)
789 write (iout,*) i-1," conformations read from DA file ",&
790 intinname(:ilen(intinname))
791 write (iout,*) jj," conformations read so far"
795 #if (defined(AIX) && !defined(JUBL))
796 call xdrfclose_(ixdrf, iret)
798 call xdrfclose(ixdrf, iret)
803 write (iout,*) "jj_old",jj_old," jj",jj
805 call write_and_send_cconf(icount,jj_old,jj,Next)
806 call MPI_Send(0,1,MPI_INTEGER,Next,570,&
807 MPI_COMM_WORLD,IERROR)
810 call write_and_send_cconf(icount,jj_old,jj,Next)
812 t_acq = tcpu() - t_acq
814 write (iout,*) "Processor",me,&
815 " time for conformation read/send",t_acq
817 ! A worker gets the confs from the master and sends them to its neighbor
819 call receive_and_pass_cconf(icount,jj_old,jj,&
821 t_acq = tcpu() - t_acq
828 write(iout,*)"A total of",ncon," conformations read."
830 allocate(enetb(1:max_ene,ncon)) !(max_ene,maxstr_proc)
832 ! Check if everyone has the same number of conformations
833 call MPI_Allgather(ncon,1,MPI_INTEGER,&
834 ntot_all(0),1,MPI_INTEGER,MPI_Comm_World,IERROR)
838 if (ncon.ne.ntot_all(i)) then
839 write (iout,*) "Number of conformations at processor",i,&
840 " differs from that at processor",me,&
848 write (iout,*) "Number of conformations read by processors"
851 write (iout,'(8i10)') i,ntot_all(i)
853 write (iout,*) "Calculation terminated."
859 1111 write(iout,*) "Error opening coordinate file ",&
860 intinname(:ilen(intinname))
863 end subroutine read_coords
864 !------------------------------------------------------------------------------
865 subroutine add_new_cconf(jjj,jj,jj_old,icount,Next)
867 use geometry_data, only: vbld,rad2deg,theta,phi,alph,omeg,deg2rad
868 use energy_data, only: itel,itype,dsc,max_ene
869 use control_data, only: symetr
870 use geometry, only: int_from_cart1
872 ! include "DIMENSIONS"
873 ! include "sizesclu.dat"
874 ! include "COMMON.CLUSTER"
875 ! include "COMMON.CONTROL"
876 ! include "COMMON.CHAIN"
877 ! include "COMMON.INTERACT"
878 ! include "COMMON.LOCAL"
879 ! include "COMMON.IOUNITS"
880 ! include "COMMON.NAMES"
881 ! include "COMMON.VAR"
882 ! include "COMMON.SBRIDGE"
883 ! include "COMMON.GEO"
884 integer :: i,j,jj,jjj,jj_old,icount,k,kk,l,ii,ib,&
885 nn,nn1,inan,Next,itj,chalen
886 real(kind=8) :: etot,energia(0:max_ene)
888 chalen=int((nct-nnt+2)/symetr)
889 call int_from_cart1(.false.)
891 if (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0) then
893 if (itel(j).ne.0 .and. itel(j-1).ne.0) then
894 write (iout,*) "Conformation",jjj,jj+1
895 write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),itel(j),&
897 write (iout,*) "The Cartesian geometry is:"
898 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
899 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
900 write (iout,*) "The internal geometry is:"
901 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
902 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
903 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
904 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
905 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
906 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
908 "This conformation WILL NOT be added to the database."
916 if (itype(j).ne.10 .and. (vbld(nres+j)-dsc(iabs(itj))) &
918 write (iout,*) "Conformation",jjj,jj+1
919 write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j)
920 write (iout,*) "The Cartesian geometry is:"
921 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
922 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
923 write (iout,*) "The internal geometry is:"
924 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
925 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
926 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
927 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
928 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
929 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
931 "This conformation WILL NOT be added to the database."
936 if (theta(j).le.0.0d0) then
938 "Zero theta angle(s) in conformation",jjj,jj+1
939 write (iout,*) "The Cartesian geometry is:"
940 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
941 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
942 write (iout,*) "The internal geometry is:"
943 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
944 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
945 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
946 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
947 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
948 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
950 "This conformation WILL NOT be added to the database."
953 if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad
957 write (iout,*) "Conformation",jjj,jj
958 write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
959 write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
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)') (vbld(k+nres),k=nnt,nct)
965 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
966 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
967 write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
968 write (iout,'(e15.5,16i5)') entfac(icount+1)
969 ! & iscore(icount+1,0)
972 call store_cconf_from_file(jj,icount)
973 if (icount.eq.maxstr_proc) then
975 write (iout,* ) "jj_old",jj_old," jj",jj
977 call write_and_send_cconf(icount,jj_old,jj,Next)
982 end subroutine add_new_cconf
983 !------------------------------------------------------------------------------
984 subroutine store_cconf_from_file(jj,icount)
986 use energy_data, only: ihpb,jhpb
988 ! include "DIMENSIONS"
989 ! include "sizesclu.dat"
990 ! include "COMMON.CLUSTER"
991 ! include "COMMON.CHAIN"
992 ! include "COMMON.SBRIDGE"
993 ! include "COMMON.INTERACT"
994 ! include "COMMON.IOUNITS"
995 ! include "COMMON.VAR"
996 integer :: i,j,jj,icount
997 ! Store the conformation that has been read in
1000 allcart(j,i,icount)=c(j,i)
1005 ihpb_all(i,icount)=ihpb(i)
1006 jhpb_all(i,icount)=jhpb(i)
1009 end subroutine store_cconf_from_file
1010 !------------------------------------------------------------------------------
1011 subroutine write_and_send_cconf(icount,jj_old,jj,Next)
1014 ! include "DIMENSIONS"
1015 ! include "sizesclu.dat"
1020 ! include "COMMON.MPI"
1022 ! include "COMMON.CHAIN"
1023 ! include "COMMON.SBRIDGE"
1024 ! include "COMMON.INTERACT"
1025 ! include "COMMON.IOUNITS"
1026 ! include "COMMON.CLUSTER"
1027 ! include "COMMON.VAR"
1028 integer :: icount,jj_old,jj,Next
1029 ! Write the structures to a scratch file
1031 ! Master sends the portion of conformations that have been read in to the neighbor
1033 write (iout,*) "Processor",me," entered WRITE_AND_SEND_CONF"
1036 call MPI_Send(icount,1,MPI_INTEGER,Next,570,MPI_COMM_WORLD,IERROR)
1037 call MPI_Send(nss_all(1),icount,MPI_INTEGER,&
1038 Next,571,MPI_COMM_WORLD,IERROR)
1039 call MPI_Send(ihpb_all(1,1),icount,MPI_INTEGER,&
1040 Next,572,MPI_COMM_WORLD,IERROR)
1041 call MPI_Send(jhpb_all(1,1),icount,MPI_INTEGER,&
1042 Next,573,MPI_COMM_WORLD,IERROR)
1043 call MPI_Send(rmstb(jj_old),icount,MPI_DOUBLE_PRECISION,&
1044 Next,577,MPI_COMM_WORLD,IERROR)
1045 call MPI_Send(entfac(jj_old),icount,MPI_DOUBLE_PRECISION,&
1046 Next,579,MPI_COMM_WORLD,IERROR)
1047 call MPI_Send(allcart(1,1,1),3*icount*2*nres,&
1048 MPI_REAL,Next,580,MPI_COMM_WORLD,IERROR)
1050 call dawrite_ccoords(jj_old,jj,icbase)
1052 end subroutine write_and_send_cconf
1053 !------------------------------------------------------------------------------
1055 subroutine receive_and_pass_cconf(icount,jj_old,jj,Previous,Next)
1059 ! include "DIMENSIONS"
1060 ! include "sizesclu.dat"
1062 integer :: IERROR !,STATUS(MPI_STATUS_SIZE)
1063 ! include "COMMON.MPI"
1064 ! include "COMMON.CHAIN"
1065 ! include "COMMON.SBRIDGE"
1066 ! include "COMMON.INTERACT"
1067 ! include "COMMON.IOUNITS"
1068 ! include "COMMON.VAR"
1069 ! include "COMMON.GEO"
1070 ! include "COMMON.CLUSTER"
1071 integer :: i,j,k,icount,jj_old,jj,Previous,Next
1074 write (iout,*) "Processor",me," entered RECEIVE_AND_PASS_CONF"
1077 do while (icount.gt.0)
1078 call MPI_Recv(icount,1,MPI_INTEGER,Previous,570,MPI_COMM_WORLD,&
1080 call MPI_Send(icount,1,MPI_INTEGER,Next,570,MPI_COMM_WORLD,&
1083 write (iout,*) "Processor",me," icount",icount
1085 if (icount.eq.0) return
1086 call MPI_Recv(nss_all(1),icount,MPI_INTEGER,&
1087 Previous,571,MPI_COMM_WORLD,STATUS,IERROR)
1088 call MPI_Send(nss_all(1),icount,MPI_INTEGER,&
1089 Next,571,MPI_COMM_WORLD,IERROR)
1090 call MPI_Recv(ihpb_all(1,1),icount,MPI_INTEGER,&
1091 Previous,572,MPI_COMM_WORLD,STATUS,IERROR)
1092 call MPI_Send(ihpb_all(1,1),icount,MPI_INTEGER,&
1093 Next,572,MPI_COMM_WORLD,IERROR)
1094 call MPI_Recv(jhpb_all(1,1),icount,MPI_INTEGER,&
1095 Previous,573,MPI_COMM_WORLD,STATUS,IERROR)
1096 call MPI_Send(jhpb_all(1,1),icount,MPI_INTEGER,&
1097 Next,573,MPI_COMM_WORLD,IERROR)
1098 call MPI_Recv(rmstb(jj_old),icount,MPI_DOUBLE_PRECISION,&
1099 Previous,577,MPI_COMM_WORLD,STATUS,IERROR)
1100 call MPI_Send(rmstb(jj_old),icount,MPI_DOUBLE_PRECISION,&
1101 Next,577,MPI_COMM_WORLD,IERROR)
1102 call MPI_Recv(entfac(jj_old),icount,MPI_DOUBLE_PRECISION,&
1103 Previous,579,MPI_COMM_WORLD,STATUS,IERROR)
1104 call MPI_Send(entfac(jj_old),icount,MPI_DOUBLE_PRECISION,&
1105 Next,579,MPI_COMM_WORLD,IERROR)
1106 call MPI_Recv(allcart(1,1,1),3*icount*2*nres,&
1107 MPI_REAL,Previous,580,MPI_COMM_WORLD,STATUS,IERROR)
1108 call MPI_Send(allcart(1,1,1),3*icount*2*nres,&
1109 MPI_REAL,Next,580,MPI_COMM_WORLD,IERROR)
1111 call dawrite_ccoords(jj_old,jj,icbase)
1114 write (iout,*) "Processor",me," received",icount," conformations"
1116 write (iout,'(8f10.4)') (allcart(l,k,i),l=1,3,k=1,nres)
1117 write (iout,'(8f10.4)')((allcart(l,k,i+nres),l=1,3,k=nnt,nct)
1118 write (iout,'(e15.5,16i5)') entfac(i)
1123 end subroutine receive_and_pass_cconf
1125 !------------------------------------------------------------------------------
1126 subroutine daread_ccoords(istart_conf,iend_conf)
1129 ! include "DIMENSIONS"
1130 ! include "sizesclu.dat"
1134 ! include "COMMON.MPI"
1136 ! include "COMMON.CHAIN"
1137 ! include "COMMON.CLUSTER"
1138 ! include "COMMON.IOUNITS"
1139 ! include "COMMON.INTERACT"
1140 ! include "COMMON.VAR"
1141 ! include "COMMON.SBRIDGE"
1142 ! include "COMMON.GEO"
1143 integer :: istart_conf,iend_conf
1144 integer :: i,j,ij,ii,iii
1146 character(len=16) :: form,acc
1147 character(len=32) :: nam
1149 ! Read conformations off a DA scratchfile.
1152 write (iout,*) "DAREAD_COORDS"
1153 write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
1154 inquire(unit=icbase,name=nam,recl=len,form=form,access=acc)
1155 write (iout,*) "len=",len," form=",form," acc=",acc
1156 write (iout,*) "nam=",nam
1159 do ii=istart_conf,iend_conf
1160 ij = ii - istart_conf + 1
1163 write (iout,*) "Reading binary file, record",iii," ii",ii
1166 read(icbase,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres),&
1167 ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
1168 nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss),&
1169 entfac(ii),rmstb(ii)
1171 write (iout,*) ii,iii,ij,entfac(ii)
1172 write (iout,'(8f10.5)') ((allcart(j,i,ij),j=1,3),i=1,nres)
1173 write (iout,'(8f10.4)') ((allcart(j,i,ij),j=1,3),&
1174 i=nnt+nres,nct+nres)
1175 write (iout,'(2e15.5)') entfac(ij)
1176 write (iout,'(16i5)') nss_all(ij),(ihpb_all(i,ij),&
1177 jhpb_all(i,ij),i=1,nss)
1182 end subroutine daread_ccoords
1183 !------------------------------------------------------------------------------
1184 subroutine dawrite_ccoords(istart_conf,iend_conf,unit_out)
1187 ! include "DIMENSIONS"
1188 ! include "sizesclu.dat"
1192 ! include "COMMON.MPI"
1194 ! include "COMMON.CHAIN"
1195 ! include "COMMON.INTERACT"
1196 ! include "COMMON.IOUNITS"
1197 ! include "COMMON.VAR"
1198 ! include "COMMON.SBRIDGE"
1199 ! include "COMMON.GEO"
1200 ! include "COMMON.CLUSTER"
1201 integer :: istart_conf,iend_conf
1202 integer :: i,j,ii,ij,iii,unit_out
1204 character(len=16) :: form,acc
1205 character(len=32) :: nam
1207 ! Write conformations to a DA scratchfile.
1210 write (iout,*) "DAWRITE_COORDS"
1211 write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
1212 write (iout,*) "lenrec",lenrec
1213 inquire(unit=unit_out,name=nam,recl=len,form=form,access=acc)
1214 write (iout,*) "len=",len," form=",form," acc=",acc
1215 write (iout,*) "nam=",nam
1218 do ii=istart_conf,iend_conf
1220 ij = ii - istart_conf + 1
1222 write (iout,*) "Writing binary file, record",iii," ii",ii
1225 write(unit_out,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres),&
1226 ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
1227 nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss_all(ij)),&
1228 entfac(ii),rmstb(ii)
1230 write (iout,'(8f10.5)') ((allcart(j,i,ij),j=1,3),i=1,nres)
1231 write (iout,'(8f10.4)') ((allcart(j,i,ij),j=1,3),i=nnt+nres,&
1233 write (iout,'(2e15.5)') entfac(ij)
1234 write (iout,'(16i5)') nss_all(ij),(ihpb(i,ij),jhpb(i,ij),i=1,&
1240 end subroutine dawrite_ccoords
1241 !-----------------------------------------------------------------------------
1243 !-----------------------------------------------------------------------------
1244 subroutine read_control
1246 ! Read molecular data
1248 use energy_data, only: rescale_mode,distchainmax,ipot !,temp0
1249 use control_data, only: titel,outpdb,outmol2,refstr,pdbref,&
1250 iscode,symetr,punch_dist,print_dist,nstart,nend,&
1251 caonly,iopt,efree,lprint_cart,lprint_int
1253 ! include 'DIMENSIONS'
1254 ! include 'sizesclu.dat'
1255 ! include 'COMMON.IOUNITS'
1256 ! include 'COMMON.TIME1'
1257 ! include 'COMMON.SBRIDGE'
1258 ! include 'COMMON.CONTROL'
1259 ! include 'COMMON.CLUSTER'
1260 ! include 'COMMON.CHAIN'
1261 ! include 'COMMON.HEADER'
1262 ! include 'COMMON.FFIELD'
1263 ! include 'COMMON.FREE'
1264 ! include 'COMMON.INTERACT'
1265 character(len=320) :: controlcard !,ucase
1267 ! include 'COMMON.INFO'
1271 read (INP,'(a80)') titel
1272 call card_concat(controlcard,.true.)
1274 call readi(controlcard,'NRES',nres,0)
1276 ! call alloc_clust_arrays
1277 allocate(rcutoff(max_cut+1)) !(max_cut+1)
1278 allocate(beta_h(maxT)) !(maxT)
1279 allocate(mult(nres)) !(maxres)
1282 call readi(controlcard,'RESCALE',rescale_mode,2)
1283 call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
1284 write (iout,*) "DISTCHAINMAX",distchainmax
1285 call readi(controlcard,'PDBOUT',outpdb,0)
1286 call readi(controlcard,'MOL2OUT',outmol2,0)
1287 refstr=(index(controlcard,'REFSTR').gt.0)
1288 write (iout,*) "REFSTR",refstr
1289 pdbref=(index(controlcard,'PDBREF').gt.0)
1290 iscode=index(controlcard,'ONE_LETTER')
1291 tree=(index(controlcard,'MAKE_TREE').gt.0)
1292 min_var=(index(controlcard,'MINVAR').gt.0)
1293 plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
1294 punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
1295 call readi(controlcard,'NCUT',ncut,1)
1296 call readi(controlcard,'SYM',symetr,1)
1297 write (iout,*) 'sym', symetr
1298 call readi(controlcard,'NSTART',nstart,0)
1299 call readi(controlcard,'NEND',nend,0)
1300 call reada(controlcard,'ECUT',ecut,10.0d0)
1301 call reada(controlcard,'PROB',prob_limit,0.99d0)
1302 write (iout,*) "Probability limit",prob_limit
1303 lgrp=(index(controlcard,'LGRP').gt.0)
1304 caonly=(index(controlcard,'CA_ONLY').gt.0)
1305 print_dist=(index(controlcard,'PRINT_DIST').gt.0)
1306 call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
1307 call readi(controlcard,'IOPT',iopt,2)
1308 lside = index(controlcard,"SIDE").gt.0
1309 efree = index(controlcard,"EFREE").gt.0
1310 call readi(controlcard,'NTEMP',nT,1)
1311 write (iout,*) "nT",nT
1312 !el call reada(controlcard,'TEMP0',temp0,300.0d0) !el
1313 call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
1314 write (iout,*) "nT",nT
1315 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
1317 beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
1319 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
1320 lprint_cart=index(controlcard,"PRINT_CART") .gt.0
1321 lprint_int=index(controlcard,"PRINT_INT") .gt.0
1324 end subroutine read_control
1325 !-----------------------------------------------------------------------------
1328 ! Read molecular data.
1330 use geometry_data, only: nsup,cref,nres0,nstart_sup,nstart_seq,dc
1331 use energy_data!, only: wsc,wscp,welec,wbond,wstrain,wtor,wtor_d,&
1332 ! wang,wscloc,wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,&
1333 ! wturn3,wturn4,wturn6,wvdwpp,weights
1334 use control_data, only: titel,nstart,nend,pdbref,refstr,iscode,&
1336 use geometry, only: chainbuild,alloc_geo_arrays
1337 use energy, only: alloc_ener_arrays
1338 use control, only: rescode,setup_var,init_int_table
1339 use conform_compar, only: contact
1341 ! include 'DIMENSIONS'
1342 ! include 'COMMON.IOUNITS'
1343 ! include 'COMMON.GEO'
1344 ! include 'COMMON.VAR'
1345 ! include 'COMMON.INTERACT'
1346 ! include 'COMMON.LOCAL'
1347 ! include 'COMMON.NAMES'
1348 ! include 'COMMON.CHAIN'
1349 ! include 'COMMON.FFIELD'
1350 ! include 'COMMON.SBRIDGE'
1351 ! include 'COMMON.HEADER'
1352 ! include 'COMMON.CONTROL'
1353 ! include 'COMMON.CONTACTS'
1354 ! include 'COMMON.TIME1'
1356 ! include 'COMMON.INFO'
1358 character(len=4) :: sequence(nres) !(maxres)
1359 character(len=800) :: weightcard
1361 real(kind=8) :: x(6*nres) !(maxvar)
1362 integer :: itype_pdb(nres) !(maxres)
1368 !el allocate(weights(n_ene))
1369 allocate(weights(max_ene))
1370 call alloc_geo_arrays
1371 call alloc_ener_arrays
1372 !-----------------------------
1373 allocate(c(3,2*nres+2)) !(3,maxres2+2) maxres2=2*maxres
1374 allocate(dc(3,0:2*nres+2)) !(3,0:maxres2)
1375 allocate(itype(nres+2)) !(maxres)
1376 allocate(itel(nres+2))
1388 !--------------------------
1389 ! Read weights of the subsequent energy terms.
1390 call card_concat(weightcard,.true.)
1391 call reada(weightcard,'WSC',wsc,1.0d0)
1392 call reada(weightcard,'WLONG',wsc,wsc)
1393 call reada(weightcard,'WSCP',wscp,1.0d0)
1394 call reada(weightcard,'WELEC',welec,1.0D0)
1395 call reada(weightcard,'WVDWPP',wvdwpp,welec)
1396 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
1397 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
1398 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
1399 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
1400 call reada(weightcard,'WTURN3',wturn3,1.0D0)
1401 call reada(weightcard,'WTURN4',wturn4,1.0D0)
1402 call reada(weightcard,'WTURN6',wturn6,1.0D0)
1403 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
1404 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
1405 call reada(weightcard,'WBOND',wbond,1.0D0)
1406 call reada(weightcard,'WTOR',wtor,1.0D0)
1407 call reada(weightcard,'WTORD',wtor_d,1.0D0)
1408 call reada(weightcard,'WANG',wang,1.0D0)
1409 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
1410 call reada(weightcard,'SCAL14',scal14,0.4D0)
1411 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
1412 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
1413 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
1414 call reada(weightcard,'TEMP0',temp0,300.0d0) !!! el
1415 if (index(weightcard,'SOFT').gt.0) ipot=6
1416 ! 12/1/95 Added weight for the multi-body term WCORR
1417 call reada(weightcard,'WCORRH',wcorr,1.0D0)
1418 if (wcorr4.gt.0.0d0) wcorr=wcorr4
1437 !el weights(19)=wsccor !!!!!!!!!!!!!!!!
1439 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
1440 wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,&
1441 wturn4,wturn6,wsccor
1442 10 format (/'Energy-term weights (unscaled):'// &
1443 'WSCC= ',f10.6,' (SC-SC)'/ &
1444 'WSCP= ',f10.6,' (SC-p)'/ &
1445 'WELEC= ',f10.6,' (p-p electr)'/ &
1446 'WVDWPP= ',f10.6,' (p-p VDW)'/ &
1447 'WBOND= ',f10.6,' (stretching)'/ &
1448 'WANG= ',f10.6,' (bending)'/ &
1449 'WSCLOC= ',f10.6,' (SC local)'/ &
1450 'WTOR= ',f10.6,' (torsional)'/ &
1451 'WTORD= ',f10.6,' (double torsional)'/ &
1452 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ &
1453 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ &
1454 'WCORR4= ',f10.6,' (multi-body 4th order)'/ &
1455 'WCORR5= ',f10.6,' (multi-body 5th order)'/ &
1456 'WCORR6= ',f10.6,' (multi-body 6th order)'/ &
1457 'WTURN3= ',f10.6,' (turns, 3rd order)'/ &
1458 'WTURN4= ',f10.6,' (turns, 4th order)'/ &
1459 'WTURN6= ',f10.6,' (turns, 6th order)'/ &
1460 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
1462 if (wcorr4.gt.0.0d0) then
1463 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',&
1464 'between contact pairs of peptide groups'
1465 write (iout,'(2(a,f5.3/))') &
1466 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,&
1467 'Range of quenching the correlation terms:',2*delt_corr
1468 else if (wcorr.gt.0.0d0) then
1469 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',&
1470 'between contact pairs of peptide groups'
1472 write (iout,'(a,f8.3)') &
1473 'Scaling factor of 1,4 SC-p interactions:',scal14
1474 write (iout,'(a,f8.3)') &
1475 'General scaling factor of SC-p interactions:',scalscp
1476 r0_corr=cutoff_corr-delt_corr
1478 aad(i,1)=scalscp*aad(i,1)
1479 aad(i,2)=scalscp*aad(i,2)
1480 bad(i,1)=scalscp*bad(i,1)
1481 bad(i,2)=scalscp*bad(i,2)
1485 print *,'indpdb=',indpdb,' pdbref=',pdbref
1487 ! Read sequence if not taken from the pdb file.
1488 if (iscode.gt.0) then
1489 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
1491 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
1493 ! Convert sequence to numeric code
1495 itype(i)=rescode(i,sequence(i),iscode)
1498 print '(20i4)',(itype(i),i=1,nres)
1502 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
1504 if (itype(i).eq.ntyp1) then
1508 else if (iabs(itype(i+1)).ne.20) then
1510 else if (iabs(itype(i)).ne.20) then
1517 write (iout,*) "ITEL"
1519 write (iout,*) i,itype(i),itel(i)
1522 print *,'Call Read_Bridge.'
1526 print *,'NNT=',NNT,' NCT=',NCT
1527 if (itype(1).eq.ntyp1) nnt=2
1528 if (itype(nres).eq.ntyp1) nct=nct-1
1529 if (nstart.lt.nnt) nstart=nnt
1530 if (nend.gt.nct .or. nend.eq.0) nend=nct
1531 write (iout,*) "nstart",nstart," nend",nend
1534 ! read(inp,'(a)') pdbfile
1535 ! write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
1536 ! open(ipdbin,file=pdbfile,status='old',err=33)
1538 ! 33 write (iout,'(a)') 'Error opening PDB file.'
1541 ! print *,'Begin reading pdb data'
1543 ! print *,'Finished reading pdb data'
1544 ! write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
1546 ! itype_pdb(i)=itype(i)
1549 ! write (iout,'(a,i3)') 'nsup=',nsup
1551 ! if (nsup.le.(nct-nnt+1)) then
1552 ! do i=0,nct-nnt+1-nsup
1553 ! if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1558 ! write (iout,'(a)')
1559 ! & 'Error - sequences to be superposed do not match.'
1562 ! do i=0,nsup-(nct-nnt+1)
1563 ! if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1565 ! nstart_sup=nstart_sup+i
1570 ! write (iout,'(a)')
1571 ! & 'Error - sequences to be superposed do not match.'
1574 ! write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1575 ! & ' nstart_seq=',nstart_seq
1577 write(iout,*)"przed ini_int_tab"
1579 write(iout,*)"po ini_int_tab"
1580 write(iout,*)"przed setup var"
1582 write(iout,*)"po setup var"
1583 write (iout,*) "molread: REFSTR",refstr
1585 if (.not.pdbref) then
1586 call read_angles(inp,*38)
1588 38 write (iout,'(a)') 'Error reading reference structure.'
1590 call mp_stopall(Error_Msg)
1592 stop 'Error reading reference structure'
1601 cref(j,i,kkk)=c(j,i)
1605 call contact(.true.,ncont_ref,icont_ref)
1608 end subroutine molread
1609 !-----------------------------------------------------------------------------
1610 subroutine openunits
1612 ! include 'DIMENSIONS'
1613 use control_data, only: from_cx,from_bx,from_cart
1617 character(len=3) :: liczba
1618 ! include "COMMON.MPI"
1620 ! include 'COMMON.IOUNITS'
1621 ! include 'COMMON.CONTROL'
1622 integer :: lenpre,lenpot !,ilen
1624 character(len=16) :: cformat,cprint
1625 ! character(len=16) ucase
1626 integer :: lenint,lenout
1627 call getenv('INPUT',prefix)
1628 call getenv('OUTPUT',prefout)
1629 call getenv('INTIN',prefintin)
1630 call getenv('COORD',cformat)
1631 call getenv('PRINTCOOR',cprint)
1632 call getenv('SCRATCHDIR',scratchdir)
1635 if (index(ucase(cformat),'CX').gt.0) then
1641 lenout=ilen(prefout)
1642 lenint=ilen(prefintin)
1643 ! Get the names and open the input files
1644 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
1646 write (liczba,'(bz,i3.3)') me
1647 outname=prefout(:lenout)//'_clust.out_'//liczba
1649 outname=prefout(:lenout)//'_clust.out'
1652 intinname=prefintin(:lenint)//'.bx'
1653 else if (from_cx) then
1654 intinname=prefintin(:lenint)//'.cx'
1656 intinname=prefintin(:lenint)//'.int'
1658 rmsname=prefintin(:lenint)//'.rms'
1659 open (jplot,file=prefout(:ilen(prefout))//'.tex',&
1661 open (jrms,file=rmsname,status='unknown')
1662 open(iout,file=outname,status='unknown')
1663 ! Get parameter filenames and open the parameter files.
1664 call getenv('BONDPAR',bondname)
1665 open (ibond,file=bondname,status='old')
1666 call getenv('THETPAR',thetname)
1667 open (ithep,file=thetname,status='old')
1668 call getenv('ROTPAR',rotname)
1669 open (irotam,file=rotname,status='old')
1670 call getenv('TORPAR',torname)
1671 open (itorp,file=torname,status='old')
1672 call getenv('TORDPAR',tordname)
1673 open (itordp,file=tordname,status='old')
1674 call getenv('FOURIER',fouriername)
1675 open (ifourier,file=fouriername,status='old')
1676 call getenv('ELEPAR',elename)
1677 open (ielep,file=elename,status='old')
1678 call getenv('SIDEPAR',sidename)
1679 open (isidep,file=sidename,status='old')
1680 call getenv('SIDEP',sidepname)
1681 open (isidep1,file=sidepname,status="old")
1682 call getenv('SCCORPAR',sccorname)
1683 open (isccor,file=sccorname,status="old")
1686 ! 8/9/01 In the newest version SCp interaction constants are read from a file
1687 ! Use -DOLDSCP to use hard-coded constants instead.
1689 call getenv('SCPPAR',scpname)
1690 open (iscpp,file=scpname,status='old')
1693 end subroutine openunits
1694 !-----------------------------------------------------------------------------
1696 !-----------------------------------------------------------------------------
1697 subroutine pdboutC(etot,rmsd,tytul)
1699 use energy_data, only: ihpb,jhpb,itype
1700 ! implicit real*8 (a-h,o-z)
1701 ! include 'DIMENSIONS'
1702 ! include 'COMMON.CONTROL'
1703 ! include 'COMMON.CHAIN'
1704 ! include 'COMMON.INTERACT'
1705 ! include 'COMMON.NAMES'
1706 ! include 'COMMON.IOUNITS'
1707 ! include 'COMMON.HEADER'
1708 ! include 'COMMON.SBRIDGE'
1709 ! include 'COMMON.TEMPFAC'
1710 character(len=50) :: tytul
1711 character(len=1) :: chainid(10)=(/'A','B','C','D','E','F',&
1713 integer :: ica(nres)
1714 real(kind=8) :: etot,rmsd
1715 integer :: iatom,ichain,ires,i,j,iti
1717 write (ipdb,'(3a,1pe15.5,a,0pf7.2)') 'REMARK ',tytul(:20),&
1718 ' ENERGY ',etot,' RMS ',rmsd
1724 if (iti.eq.ntyp1) then
1727 write (ipdb,'(a)') 'TER'
1732 write (ipdb,10) iatom,restyp(iti),chainid(ichain),&
1733 ires,(c(j,i),j=1,3),1.0d0,tempfac(1,i)
1736 write (ipdb,20) iatom,restyp(iti),chainid(ichain),&
1737 ires,(c(j,nres+i),j=1,3),1.0d0,tempfac(2,i)
1741 write (ipdb,'(a)') 'TER'
1743 if (itype(i).eq.ntyp1) cycle
1744 if (itype(i).eq.10 .and. itype(i+1).ne.ntyp1) then
1745 write (ipdb,30) ica(i),ica(i+1)
1746 else if (itype(i).ne.10 .and. itype(i+1).ne.ntyp1) then
1747 write (ipdb,30) ica(i),ica(i+1),ica(i)+1
1748 else if (itype(i).ne.10 .and. itype(i+1).eq.ntyp1) then
1749 write (ipdb,30) ica(i),ica(i)+1
1752 if (itype(nct).ne.10) then
1753 write (ipdb,30) ica(nct),ica(nct)+1
1756 write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
1758 write (ipdb,'(a6)') 'ENDMDL'
1759 10 FORMAT ('ATOM',I7,' CA ',A3,1X,A1,I4,4X,3F8.3,2f6.2)
1760 20 FORMAT ('ATOM',I7,' CB ',A3,1X,A1,I4,4X,3F8.3,2f6.2)
1761 30 FORMAT ('CONECT',8I5)
1763 end subroutine pdboutC
1764 !-----------------------------------------------------------------------------
1765 subroutine cartout(igr,i,etot,free,rmsd,plik)
1766 ! implicit real*8 (a-h,o-z)
1767 ! include 'DIMENSIONS'
1768 ! include 'sizesclu.dat'
1769 ! include 'COMMON.IOUNITS'
1770 ! include 'COMMON.CHAIN'
1771 ! include 'COMMON.VAR'
1772 ! include 'COMMON.LOCAL'
1773 ! include 'COMMON.INTERACT'
1774 ! include 'COMMON.NAMES'
1775 ! include 'COMMON.GEO'
1776 ! include 'COMMON.CLUSTER'
1777 integer :: igr,i,j,k
1778 real(kind=8) :: etot,free,rmsd
1779 character(len=80) :: plik
1780 open (igeom,file=plik,position='append')
1781 write (igeom,'(2e15.5,f10.5,$)') etot,free,rmsd
1782 write (igeom,'(i4,$)') &
1783 nss_all(i),(ihpb_all(j,i),jhpb_all(j,i),j=1,nss_all(i))
1784 write (igeom,'(i10)') iscore(i)
1785 write (igeom,'(8f10.5)') &
1786 ((allcart(k,j,i),k=1,3),j=1,nres),&
1787 ((allcart(k,j+nres,i),k=1,3),j=nnt,nct)
1789 end subroutine cartout
1790 !------------------------------------------------------------------------------
1791 ! subroutine alloc_clust_arrays(n_conf)
1796 ! allocate(diss(maxdist)) !(maxdist)
1797 !el allocate(energy(0:maxconf),totfree(0:maxconf)) !(0:maxconf)
1798 ! allocatable :: enetb !(max_ene,maxstr_proc)
1799 !el allocate(entfac(maxconf)) !(maxconf)
1800 ! allocatable :: totfree_gr !(maxgr)
1801 !el allocate(rcutoff(max_cut+1)) !(max_cut+1)
1803 ! allocatable :: licz,iass !(maxgr)
1804 ! allocatable :: nconf !(maxgr,maxingr)
1805 ! allocatable :: iass_tot !(maxgr,max_cut)
1806 ! allocatable :: list_conf !(maxconf)
1808 !el allocatable :: allcart !(3,maxres2,maxstr_proc)
1809 !el allocate(rmstb(maxconf)) !(maxconf)
1810 !el allocate(mult(nres)) !(maxres)
1811 !el allocatable :: nss_all !(maxstr_proc)
1812 !el allocatable :: ihpb_all,jhpb_all !(maxss,maxstr_proc)
1813 ! allocate(icc(n_conf),iscore(n_conf)) !(maxconf)
1816 ! allocatable :: tempfac !(2,maxres)
1819 !el allocate(beta_h(maxT)) !(maxT)
1821 ! end subroutine alloc_clust_arrays
1822 !-----------------------------------------------------------------------------
1823 !-----------------------------------------------------------------------------