2 !-----------------------------------------------------------------------------
6 use io_base !, only: ilen
7 use geometry_data, only: nres,c,nres_molec,boxxsize,boxysize,boxzsize
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,molnum
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,mnum
886 real(kind=8) :: etot,energia(0:max_ene)
888 chalen=int((nct-nnt+2)/symetr)
889 call int_from_cart1(.false.)
893 if (itype(j,mnum).eq.ntyp1_molec(mnum)) cycle
894 if (itype(j+1,molnum(j+1)).eq.ntyp1_molec(molnum(j+1))) cycle
896 if ((vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0).and.(mnum.ne.5)) then
898 if (itel(j).ne.0 .and. itel(j-1).ne.0) then
899 write (iout,*) "Conformation",jjj,jj+1
900 write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),itel(j),&
902 write (iout,*) "The Cartesian geometry is:"
903 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
904 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
905 write (iout,*) "The internal geometry is:"
906 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
907 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
908 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
909 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
910 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
911 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
913 "This conformation WILL NOT be added to the database."
923 if (itype(j,1).ne.10 .and. (itype(j,mnum).ne.ntyp1_molec(mnum)) &
924 .and. (vbld(nres+j)-dsc(iabs(itj))) &
926 write (iout,*) "Conformation",jjj,jj+1
927 write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j)
928 write (iout,*) "The Cartesian geometry is:"
929 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
930 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
931 write (iout,*) "The internal geometry is:"
932 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
933 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
934 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
935 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
936 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
937 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
939 "This conformation WILL NOT be added to the database."
944 if (theta(j).le.0.0d0) then
946 "Zero theta angle(s) in conformation",jjj,jj+1
947 write (iout,*) "The Cartesian geometry is:"
948 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
949 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
950 write (iout,*) "The internal geometry is:"
951 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
952 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
953 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
954 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
955 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
956 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
958 "This conformation WILL NOT be added to the database."
961 if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad
965 write (iout,*) "Conformation",jjj,jj
966 write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
967 write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
968 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
969 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
970 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
971 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
972 write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct)
973 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
974 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
975 write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
976 write (iout,'(e15.5,16i5)') entfac(icount+1)
977 ! & iscore(icount+1,0)
980 call store_cconf_from_file(jj,icount)
981 if (icount.eq.maxstr_proc) then
983 write (iout,* ) "jj_old",jj_old," jj",jj
985 call write_and_send_cconf(icount,jj_old,jj,Next)
990 end subroutine add_new_cconf
991 !------------------------------------------------------------------------------
992 subroutine store_cconf_from_file(jj,icount)
994 use energy_data, only: ihpb,jhpb
996 ! include "DIMENSIONS"
997 ! include "sizesclu.dat"
998 ! include "COMMON.CLUSTER"
999 ! include "COMMON.CHAIN"
1000 ! include "COMMON.SBRIDGE"
1001 ! include "COMMON.INTERACT"
1002 ! include "COMMON.IOUNITS"
1003 ! include "COMMON.VAR"
1004 integer :: i,j,jj,icount
1005 ! Store the conformation that has been read in
1008 allcart(j,i,icount)=c(j,i)
1013 ihpb_all(i,icount)=ihpb(i)
1014 jhpb_all(i,icount)=jhpb(i)
1017 end subroutine store_cconf_from_file
1018 !------------------------------------------------------------------------------
1019 subroutine write_and_send_cconf(icount,jj_old,jj,Next)
1022 ! include "DIMENSIONS"
1023 ! include "sizesclu.dat"
1028 ! include "COMMON.MPI"
1030 ! include "COMMON.CHAIN"
1031 ! include "COMMON.SBRIDGE"
1032 ! include "COMMON.INTERACT"
1033 ! include "COMMON.IOUNITS"
1034 ! include "COMMON.CLUSTER"
1035 ! include "COMMON.VAR"
1036 integer :: icount,jj_old,jj,Next
1037 ! Write the structures to a scratch file
1039 ! Master sends the portion of conformations that have been read in to the neighbor
1041 write (iout,*) "Processor",me," entered WRITE_AND_SEND_CONF"
1044 call MPI_Send(icount,1,MPI_INTEGER,Next,570,MPI_COMM_WORLD,IERROR)
1045 call MPI_Send(nss_all(1),icount,MPI_INTEGER,&
1046 Next,571,MPI_COMM_WORLD,IERROR)
1047 call MPI_Send(ihpb_all(1,1),icount,MPI_INTEGER,&
1048 Next,572,MPI_COMM_WORLD,IERROR)
1049 call MPI_Send(jhpb_all(1,1),icount,MPI_INTEGER,&
1050 Next,573,MPI_COMM_WORLD,IERROR)
1051 call MPI_Send(rmstb(jj_old),icount,MPI_DOUBLE_PRECISION,&
1052 Next,577,MPI_COMM_WORLD,IERROR)
1053 call MPI_Send(entfac(jj_old),icount,MPI_DOUBLE_PRECISION,&
1054 Next,579,MPI_COMM_WORLD,IERROR)
1055 call MPI_Send(allcart(1,1,1),3*icount*2*nres,&
1056 MPI_REAL,Next,580,MPI_COMM_WORLD,IERROR)
1058 call dawrite_ccoords(jj_old,jj,icbase)
1060 end subroutine write_and_send_cconf
1061 !------------------------------------------------------------------------------
1063 subroutine receive_and_pass_cconf(icount,jj_old,jj,Previous,Next)
1067 ! include "DIMENSIONS"
1068 ! include "sizesclu.dat"
1070 integer :: IERROR !,STATUS(MPI_STATUS_SIZE)
1071 ! include "COMMON.MPI"
1072 ! include "COMMON.CHAIN"
1073 ! include "COMMON.SBRIDGE"
1074 ! include "COMMON.INTERACT"
1075 ! include "COMMON.IOUNITS"
1076 ! include "COMMON.VAR"
1077 ! include "COMMON.GEO"
1078 ! include "COMMON.CLUSTER"
1079 integer :: i,j,k,icount,jj_old,jj,Previous,Next
1082 write (iout,*) "Processor",me," entered RECEIVE_AND_PASS_CONF"
1085 do while (icount.gt.0)
1086 call MPI_Recv(icount,1,MPI_INTEGER,Previous,570,MPI_COMM_WORLD,&
1088 call MPI_Send(icount,1,MPI_INTEGER,Next,570,MPI_COMM_WORLD,&
1091 write (iout,*) "Processor",me," icount",icount
1093 if (icount.eq.0) return
1094 call MPI_Recv(nss_all(1),icount,MPI_INTEGER,&
1095 Previous,571,MPI_COMM_WORLD,STATUS,IERROR)
1096 call MPI_Send(nss_all(1),icount,MPI_INTEGER,&
1097 Next,571,MPI_COMM_WORLD,IERROR)
1098 call MPI_Recv(ihpb_all(1,1),icount,MPI_INTEGER,&
1099 Previous,572,MPI_COMM_WORLD,STATUS,IERROR)
1100 call MPI_Send(ihpb_all(1,1),icount,MPI_INTEGER,&
1101 Next,572,MPI_COMM_WORLD,IERROR)
1102 call MPI_Recv(jhpb_all(1,1),icount,MPI_INTEGER,&
1103 Previous,573,MPI_COMM_WORLD,STATUS,IERROR)
1104 call MPI_Send(jhpb_all(1,1),icount,MPI_INTEGER,&
1105 Next,573,MPI_COMM_WORLD,IERROR)
1106 call MPI_Recv(rmstb(jj_old),icount,MPI_DOUBLE_PRECISION,&
1107 Previous,577,MPI_COMM_WORLD,STATUS,IERROR)
1108 call MPI_Send(rmstb(jj_old),icount,MPI_DOUBLE_PRECISION,&
1109 Next,577,MPI_COMM_WORLD,IERROR)
1110 call MPI_Recv(entfac(jj_old),icount,MPI_DOUBLE_PRECISION,&
1111 Previous,579,MPI_COMM_WORLD,STATUS,IERROR)
1112 call MPI_Send(entfac(jj_old),icount,MPI_DOUBLE_PRECISION,&
1113 Next,579,MPI_COMM_WORLD,IERROR)
1114 call MPI_Recv(allcart(1,1,1),3*icount*2*nres,&
1115 MPI_REAL,Previous,580,MPI_COMM_WORLD,STATUS,IERROR)
1116 call MPI_Send(allcart(1,1,1),3*icount*2*nres,&
1117 MPI_REAL,Next,580,MPI_COMM_WORLD,IERROR)
1119 call dawrite_ccoords(jj_old,jj,icbase)
1122 write (iout,*) "Processor",me," received",icount," conformations"
1124 write (iout,'(8f10.4)') (allcart(l,k,i),l=1,3,k=1,nres)
1125 write (iout,'(8f10.4)')((allcart(l,k,i+nres),l=1,3,k=nnt,nct)
1126 write (iout,'(e15.5,16i5)') entfac(i)
1131 end subroutine receive_and_pass_cconf
1133 !------------------------------------------------------------------------------
1134 subroutine daread_ccoords(istart_conf,iend_conf)
1137 ! include "DIMENSIONS"
1138 ! include "sizesclu.dat"
1142 ! include "COMMON.MPI"
1144 ! include "COMMON.CHAIN"
1145 ! include "COMMON.CLUSTER"
1146 ! include "COMMON.IOUNITS"
1147 ! include "COMMON.INTERACT"
1148 ! include "COMMON.VAR"
1149 ! include "COMMON.SBRIDGE"
1150 ! include "COMMON.GEO"
1151 integer :: istart_conf,iend_conf
1152 integer :: i,j,ij,ii,iii
1154 character(len=16) :: form,acc
1155 character(len=32) :: nam
1157 ! Read conformations off a DA scratchfile.
1160 write (iout,*) "DAREAD_COORDS"
1161 write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
1162 inquire(unit=icbase,name=nam,recl=len,form=form,access=acc)
1163 write (iout,*) "len=",len," form=",form," acc=",acc
1164 write (iout,*) "nam=",nam
1167 do ii=istart_conf,iend_conf
1168 ij = ii - istart_conf + 1
1171 write (iout,*) "Reading binary file, record",iii," ii",ii
1174 read(icbase,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres),&
1175 ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
1176 nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss),&
1177 entfac(ii),rmstb(ii)
1179 write (iout,*) ii,iii,ij,entfac(ii)
1180 write (iout,'(8f10.5)') ((allcart(j,i,ij),j=1,3),i=1,nres)
1181 write (iout,'(8f10.4)') ((allcart(j,i,ij),j=1,3),&
1182 i=nnt+nres,nct+nres)
1183 write (iout,'(2e15.5)') entfac(ij)
1184 write (iout,'(16i5)') nss_all(ij),(ihpb_all(i,ij),&
1185 jhpb_all(i,ij),i=1,nss)
1190 end subroutine daread_ccoords
1191 !------------------------------------------------------------------------------
1192 subroutine dawrite_ccoords(istart_conf,iend_conf,unit_out)
1195 ! include "DIMENSIONS"
1196 ! include "sizesclu.dat"
1200 ! include "COMMON.MPI"
1202 ! include "COMMON.CHAIN"
1203 ! include "COMMON.INTERACT"
1204 ! include "COMMON.IOUNITS"
1205 ! include "COMMON.VAR"
1206 ! include "COMMON.SBRIDGE"
1207 ! include "COMMON.GEO"
1208 ! include "COMMON.CLUSTER"
1209 integer :: istart_conf,iend_conf
1210 integer :: i,j,ii,ij,iii,unit_out
1212 character(len=16) :: form,acc
1213 character(len=32) :: nam
1215 ! Write conformations to a DA scratchfile.
1218 write (iout,*) "DAWRITE_COORDS"
1219 write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
1220 write (iout,*) "lenrec",lenrec
1221 inquire(unit=unit_out,name=nam,recl=len,form=form,access=acc)
1222 write (iout,*) "len=",len," form=",form," acc=",acc
1223 write (iout,*) "nam=",nam
1226 do ii=istart_conf,iend_conf
1228 ij = ii - istart_conf + 1
1230 write (iout,*) "Writing binary file, record",iii," ii",ii
1233 write(unit_out,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres),&
1234 ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
1235 nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss_all(ij)),&
1236 entfac(ii),rmstb(ii)
1238 write (iout,'(8f10.5)') ((allcart(j,i,ij),j=1,3),i=1,nres)
1239 write (iout,'(8f10.4)') ((allcart(j,i,ij),j=1,3),i=nnt+nres,&
1241 write (iout,'(2e15.5)') entfac(ij)
1242 write (iout,'(16i5)') nss_all(ij),(ihpb(i,ij),jhpb(i,ij),i=1,&
1248 end subroutine dawrite_ccoords
1249 !-----------------------------------------------------------------------------
1251 !-----------------------------------------------------------------------------
1252 subroutine read_control
1254 ! Read molecular data
1256 use energy_data, only: rescale_mode,distchainmax,ipot !,temp0
1257 use control_data, only: titel,outpdb,outmol2,refstr,pdbref,&
1258 iscode,symetr,punch_dist,print_dist,nstart,nend,&
1259 caonly,iopt,efree,lprint_cart,lprint_int,rlamb_ele,&
1262 ! include 'DIMENSIONS'
1263 ! include 'sizesclu.dat'
1264 ! include 'COMMON.IOUNITS'
1265 ! include 'COMMON.TIME1'
1266 ! include 'COMMON.SBRIDGE'
1267 ! include 'COMMON.CONTROL'
1268 ! include 'COMMON.CLUSTER'
1269 ! include 'COMMON.CHAIN'
1270 ! include 'COMMON.HEADER'
1271 ! include 'COMMON.FFIELD'
1272 ! include 'COMMON.FREE'
1273 ! include 'COMMON.INTERACT'
1274 character(len=320) :: controlcard !,ucase
1276 ! include 'COMMON.INFO'
1280 read (INP,'(a80)') titel
1281 call card_concat(controlcard,.true.)
1283 call readi(controlcard,'NRES',nres_molec(1),0)
1284 call readi(controlcard,"NRES_NUCL",nres_molec(2),0)
1285 call readi(controlcard,"NRES_CAT",nres_molec(5),0)
1288 nres=nres_molec(i)+nres
1291 ! call alloc_clust_arrays
1292 allocate(rcutoff(max_cut+1)) !(max_cut+1)
1293 allocate(beta_h(maxT)) !(maxT)
1294 allocate(mult(nres)) !(maxres)
1297 call readi(controlcard,'RESCALE',rescale_mode,2)
1298 call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
1299 write (iout,*) "DISTCHAINMAX",distchainmax
1300 call readi(controlcard,'PDBOUT',outpdb,0)
1301 call readi(controlcard,'MOL2OUT',outmol2,0)
1302 refstr=(index(controlcard,'REFSTR').gt.0)
1303 write (iout,*) "REFSTR",refstr
1304 pdbref=(index(controlcard,'PDBREF').gt.0)
1305 iscode=index(controlcard,'ONE_LETTER')
1306 tree=(index(controlcard,'MAKE_TREE').gt.0)
1307 min_var=(index(controlcard,'MINVAR').gt.0)
1308 plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
1309 punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
1310 call readi(controlcard,'NCUT',ncut,1)
1311 call readi(controlcard,'SYM',symetr,1)
1312 write (iout,*) 'sym', symetr
1313 call readi(controlcard,'NSTART',nstart,0)
1314 call reada(controlcard,'BOXX',boxxsize,100.0d0)
1315 call reada(controlcard,'BOXY',boxysize,100.0d0)
1316 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
1317 ! ions=index(controlcard,"IONS").gt.0
1318 call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
1319 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
1320 write(iout,*) "R_CUT_ELE=",r_cut_ele
1321 call readi(controlcard,'NEND',nend,0)
1322 call reada(controlcard,'ECUT',ecut,10.0d0)
1323 call reada(controlcard,'PROB',prob_limit,0.99d0)
1324 write (iout,*) "Probability limit",prob_limit
1325 lgrp=(index(controlcard,'LGRP').gt.0)
1326 caonly=(index(controlcard,'CA_ONLY').gt.0)
1327 print_dist=(index(controlcard,'PRINT_DIST').gt.0)
1328 call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
1329 call readi(controlcard,'IOPT',iopt,2)
1330 lside = index(controlcard,"SIDE").gt.0
1331 efree = index(controlcard,"EFREE").gt.0
1332 call readi(controlcard,'NTEMP',nT,1)
1333 write (iout,*) "nT",nT
1334 !el call reada(controlcard,'TEMP0',temp0,300.0d0) !el
1335 call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
1336 write (iout,*) "nT",nT
1337 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
1339 beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
1341 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
1342 lprint_cart=index(controlcard,"PRINT_CART") .gt.0
1343 lprint_int=index(controlcard,"PRINT_INT") .gt.0
1346 end subroutine read_control
1347 !-----------------------------------------------------------------------------
1350 ! Read molecular data.
1352 use geometry_data, only: nsup,cref,nres0,nstart_sup,nstart_seq,dc
1353 use energy_data!, only: wsc,wscp,welec,wbond,wstrain,wtor,wtor_d,&
1354 ! wang,wscloc,wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,&
1355 ! wturn3,wturn4,wturn6,wvdwpp,weights
1356 use control_data, only: titel,nstart,nend,pdbref,refstr,iscode,&
1358 use geometry, only: chainbuild,alloc_geo_arrays
1359 use energy, only: alloc_ener_arrays
1360 use control, only: rescode,setup_var,init_int_table
1361 use conform_compar, only: contact
1363 ! include 'DIMENSIONS'
1364 ! include 'COMMON.IOUNITS'
1365 ! include 'COMMON.GEO'
1366 ! include 'COMMON.VAR'
1367 ! include 'COMMON.INTERACT'
1368 ! include 'COMMON.LOCAL'
1369 ! include 'COMMON.NAMES'
1370 ! include 'COMMON.CHAIN'
1371 ! include 'COMMON.FFIELD'
1372 ! include 'COMMON.SBRIDGE'
1373 ! include 'COMMON.HEADER'
1374 ! include 'COMMON.CONTROL'
1375 ! include 'COMMON.CONTACTS'
1376 ! include 'COMMON.TIME1'
1378 ! include 'COMMON.INFO'
1380 character(len=4) :: sequence(nres) !(maxres)
1381 character(len=800) :: weightcard
1383 real(kind=8) :: x(6*nres) !(maxvar)
1384 integer :: itype_pdb(nres) !(maxres)
1386 integer :: i,j,kkk,mnum
1390 !el allocate(weights(n_ene))
1391 allocate(weights(max_ene))
1392 call alloc_geo_arrays
1393 call alloc_ener_arrays
1394 !-----------------------------
1395 allocate(c(3,2*nres+2)) !(3,maxres2+2) maxres2=2*maxres
1396 allocate(dc(3,0:2*nres+2)) !(3,0:maxres2)
1397 allocate(itype(nres+2,5)) !(maxres)
1398 allocate(molnum(nres+1))
1399 allocate(itel(nres+2))
1413 !--------------------------
1414 ! Read weights of the subsequent energy terms.
1415 call card_concat(weightcard,.true.)
1416 call reada(weightcard,'WSC',wsc,1.0d0)
1417 call reada(weightcard,'WLONG',wsc,wsc)
1418 call reada(weightcard,'WSCP',wscp,1.0d0)
1419 call reada(weightcard,'WELEC',welec,1.0D0)
1420 call reada(weightcard,'WVDWPP',wvdwpp,welec)
1421 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
1422 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
1423 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
1424 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
1425 call reada(weightcard,'WTURN3',wturn3,1.0D0)
1426 call reada(weightcard,'WTURN4',wturn4,1.0D0)
1427 call reada(weightcard,'WTURN6',wturn6,1.0D0)
1428 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
1429 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
1430 call reada(weightcard,'WBOND',wbond,1.0D0)
1431 call reada(weightcard,'WTOR',wtor,1.0D0)
1432 call reada(weightcard,'WTORD',wtor_d,1.0D0)
1433 call reada(weightcard,'WANG',wang,1.0D0)
1434 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
1435 call reada(weightcard,'SCAL14',scal14,0.4D0)
1436 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
1437 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
1438 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
1439 call reada(weightcard,'WCATCAT',wcatcat,0.0d0)
1440 call reada(weightcard,'WCATPROT',wcatprot,0.0d0)
1441 call reada(weightcard,'TEMP0',temp0,300.0d0) !!! el
1442 if (index(weightcard,'SOFT').gt.0) ipot=6
1443 ! 12/1/95 Added weight for the multi-body term WCORR
1444 call reada(weightcard,'WCORRH',wcorr,1.0D0)
1445 if (wcorr4.gt.0.0d0) wcorr=wcorr4
1464 !el weights(19)=wsccor !!!!!!!!!!!!!!!!
1466 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
1467 wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,&
1468 wturn4,wturn6,wsccor
1469 10 format (/'Energy-term weights (unscaled):'// &
1470 'WSCC= ',f10.6,' (SC-SC)'/ &
1471 'WSCP= ',f10.6,' (SC-p)'/ &
1472 'WELEC= ',f10.6,' (p-p electr)'/ &
1473 'WVDWPP= ',f10.6,' (p-p VDW)'/ &
1474 'WBOND= ',f10.6,' (stretching)'/ &
1475 'WANG= ',f10.6,' (bending)'/ &
1476 'WSCLOC= ',f10.6,' (SC local)'/ &
1477 'WTOR= ',f10.6,' (torsional)'/ &
1478 'WTORD= ',f10.6,' (double torsional)'/ &
1479 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ &
1480 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ &
1481 'WCORR4= ',f10.6,' (multi-body 4th order)'/ &
1482 'WCORR5= ',f10.6,' (multi-body 5th order)'/ &
1483 'WCORR6= ',f10.6,' (multi-body 6th order)'/ &
1484 'WTURN3= ',f10.6,' (turns, 3rd order)'/ &
1485 'WTURN4= ',f10.6,' (turns, 4th order)'/ &
1486 'WTURN6= ',f10.6,' (turns, 6th order)'/ &
1487 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
1489 if (wcorr4.gt.0.0d0) then
1490 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',&
1491 'between contact pairs of peptide groups'
1492 write (iout,'(2(a,f5.3/))') &
1493 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,&
1494 'Range of quenching the correlation terms:',2*delt_corr
1495 else if (wcorr.gt.0.0d0) then
1496 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',&
1497 'between contact pairs of peptide groups'
1499 write (iout,'(a,f8.3)') &
1500 'Scaling factor of 1,4 SC-p interactions:',scal14
1501 write (iout,'(a,f8.3)') &
1502 'General scaling factor of SC-p interactions:',scalscp
1503 r0_corr=cutoff_corr-delt_corr
1505 aad(i,1)=scalscp*aad(i,1)
1506 aad(i,2)=scalscp*aad(i,2)
1507 bad(i,1)=scalscp*bad(i,1)
1508 bad(i,2)=scalscp*bad(i,2)
1512 print *,'indpdb=',indpdb,' pdbref=',pdbref
1514 ! Read sequence if not taken from the pdb file.
1515 if (iscode.gt.0) then
1516 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
1518 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
1520 ! Convert sequence to numeric code
1521 do i=1,nres_molec(1)
1524 itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
1526 do i=nres_molec(1)+1,nres_molec(1)+nres_molec(2)
1529 itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
1531 do i=nres_molec(1)+nres_molec(2)+1,nres_molec(1)+nres_molec(2)+nres_molec(5)
1534 itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
1537 print '(20i4)',(itype(i,molnum(i)),i=1,nres)
1541 if (itype(i,molnum(i)).eq.ntyp1_molec(molnum(i))&
1542 .or. itype(i+1,molnum(i+1)).eq.ntyp1_molec(molnum(i+1))) then
1544 if (itype(i,molnum(i)).eq.ntyp1_molec(molnum(i))) then
1548 else if (iabs(itype(i+1,1)).ne.20) then
1550 else if (iabs(itype(i,1)).ne.20) then
1557 write (iout,*) "ITEL"
1559 write (iout,*) i,itype(i,molnum(i)),itel(i)
1562 print *,'Call Read_Bridge.'
1566 print *,'NNT=',NNT,' NCT=',NCT
1567 if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2
1568 if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1
1569 if (nstart.lt.nnt) nstart=nnt
1570 if (nend.gt.nct .or. nend.eq.0) nend=nct
1571 write (iout,*) "nstart",nstart," nend",nend
1574 ! read(inp,'(a)') pdbfile
1575 ! write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
1576 ! open(ipdbin,file=pdbfile,status='old',err=33)
1578 ! 33 write (iout,'(a)') 'Error opening PDB file.'
1581 ! print *,'Begin reading pdb data'
1583 ! print *,'Finished reading pdb data'
1584 ! write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
1586 ! itype_pdb(i)=itype(i)
1589 ! write (iout,'(a,i3)') 'nsup=',nsup
1591 ! if (nsup.le.(nct-nnt+1)) then
1592 ! do i=0,nct-nnt+1-nsup
1593 ! if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1598 ! write (iout,'(a)')
1599 ! & 'Error - sequences to be superposed do not match.'
1602 ! do i=0,nsup-(nct-nnt+1)
1603 ! if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1605 ! nstart_sup=nstart_sup+i
1610 ! write (iout,'(a)')
1611 ! & 'Error - sequences to be superposed do not match.'
1614 ! write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1615 ! & ' nstart_seq=',nstart_seq
1617 write(iout,*)"przed ini_int_tab"
1619 write(iout,*)"po ini_int_tab"
1620 write(iout,*)"przed setup var"
1622 write(iout,*)"po setup var"
1623 write (iout,*) "molread: REFSTR",refstr
1625 if (.not.pdbref) then
1626 call read_angles(inp,*38)
1628 38 write (iout,'(a)') 'Error reading reference structure.'
1630 call mp_stopall(Error_Msg)
1632 stop 'Error reading reference structure'
1641 cref(j,i,kkk)=c(j,i)
1645 call contact(.true.,ncont_ref,icont_ref)
1648 end subroutine molread
1649 !-----------------------------------------------------------------------------
1650 subroutine openunits
1652 ! include 'DIMENSIONS'
1653 use control_data, only: from_cx,from_bx,from_cart
1657 character(len=3) :: liczba
1658 ! include "COMMON.MPI"
1660 ! include 'COMMON.IOUNITS'
1661 ! include 'COMMON.CONTROL'
1662 integer :: lenpre,lenpot !,ilen
1664 character(len=16) :: cformat,cprint
1665 ! character(len=16) ucase
1666 integer :: lenint,lenout
1667 call getenv('INPUT',prefix)
1668 call getenv('OUTPUT',prefout)
1669 call getenv('INTIN',prefintin)
1670 call getenv('COORD',cformat)
1671 call getenv('PRINTCOOR',cprint)
1672 call getenv('SCRATCHDIR',scratchdir)
1675 if (index(ucase(cformat),'CX').gt.0) then
1681 lenout=ilen(prefout)
1682 lenint=ilen(prefintin)
1683 ! Get the names and open the input files
1684 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
1686 write (liczba,'(bz,i3.3)') me
1687 outname=prefout(:lenout)//'_clust.out_'//liczba
1689 outname=prefout(:lenout)//'_clust.out'
1692 intinname=prefintin(:lenint)//'.bx'
1693 else if (from_cx) then
1694 intinname=prefintin(:lenint)//'.cx'
1696 intinname=prefintin(:lenint)//'.int'
1698 rmsname=prefintin(:lenint)//'.rms'
1699 open (jplot,file=prefout(:ilen(prefout))//'.tex',&
1701 open (jrms,file=rmsname,status='unknown')
1702 open(iout,file=outname,status='unknown')
1703 ! Get parameter filenames and open the parameter files.
1704 call getenv('BONDPAR',bondname)
1705 open (ibond,file=bondname,status='old')
1706 call getenv('THETPAR',thetname)
1707 open (ithep,file=thetname,status='old')
1708 call getenv('ROTPAR',rotname)
1709 open (irotam,file=rotname,status='old')
1710 call getenv('TORPAR',torname)
1711 open (itorp,file=torname,status='old')
1712 call getenv('TORDPAR',tordname)
1713 open (itordp,file=tordname,status='old')
1714 call getenv('FOURIER',fouriername)
1715 open (ifourier,file=fouriername,status='old')
1716 call getenv('ELEPAR',elename)
1717 open (ielep,file=elename,status='old')
1718 call getenv('SIDEPAR',sidename)
1719 open (isidep,file=sidename,status='old')
1720 call getenv('SIDEP',sidepname)
1721 open (isidep1,file=sidepname,status="old")
1722 call getenv('SCCORPAR',sccorname)
1723 open (isccor,file=sccorname,status="old")
1724 call getenv('THETPAR_NUCL',thetname_nucl)
1725 open (ithep_nucl,file=thetname_nucl,status='old')
1726 call getenv('ROTPAR_NUCL',rotname_nucl)
1727 open (irotam_nucl,file=rotname_nucl,status='old')
1728 call getenv('TORPAR_NUCL',torname_nucl)
1729 open (itorp_nucl,file=torname_nucl,status='old')
1730 call getenv('TORDPAR_NUCL',tordname_nucl)
1731 open (itordp_nucl,file=tordname_nucl,status='old')
1732 call getenv('SIDEPAR_NUCL',sidename_nucl)
1733 open (isidep_nucl,file=sidename_nucl,status='old')
1734 call getenv('SIDEPAR_SCBASE',sidename_scbase)
1735 open (isidep_scbase,file=sidename_scbase,status='old')
1736 call getenv('PEPPAR_PEPBASE',pepname_pepbase)
1737 open (isidep_pepbase,file=pepname_pepbase,status='old')
1738 call getenv('SCPAR_PHOSPH',pepname_scpho)
1739 open (isidep_scpho,file=pepname_scpho,status='old')
1740 call getenv('PEPPAR_PHOSPH',pepname_peppho)
1741 open (isidep_peppho,file=pepname_peppho,status='old')
1744 call getenv('LIPTRANPAR',liptranname)
1745 open (iliptranpar,file=liptranname,status='old')
1746 call getenv('TUBEPAR',tubename)
1747 open (itube,file=tubename,status='old')
1748 call getenv('IONPAR',ionname)
1749 open (iion,file=ionname,status='old')
1753 ! 8/9/01 In the newest version SCp interaction constants are read from a file
1754 ! Use -DOLDSCP to use hard-coded constants instead.
1756 call getenv('SCPPAR',scpname)
1757 open (iscpp,file=scpname,status='old')
1760 end subroutine openunits
1761 !-----------------------------------------------------------------------------
1763 !-----------------------------------------------------------------------------
1764 subroutine pdboutC(etot,rmsd,tytul)
1766 use energy_data, only: ihpb,jhpb,itype,molnum
1767 ! implicit real*8 (a-h,o-z)
1768 ! include 'DIMENSIONS'
1769 ! include 'COMMON.CONTROL'
1770 ! include 'COMMON.CHAIN'
1771 ! include 'COMMON.INTERACT'
1772 ! include 'COMMON.NAMES'
1773 ! include 'COMMON.IOUNITS'
1774 ! include 'COMMON.HEADER'
1775 ! include 'COMMON.SBRIDGE'
1776 ! include 'COMMON.TEMPFAC'
1777 character(len=50) :: tytul
1778 character(len=1) :: chainid(10)=(/'A','B','C','D','E','F',&
1780 integer :: ica(nres)
1781 real(kind=8) :: etot,rmsd
1782 integer :: iatom,ichain,ires,i,j,iti,mnum,mnum1,boxxxshift(3)
1783 real(kind=8) :: boxxxx(3)
1784 write (ipdb,'(3a,1pe15.5,a,0pf7.2)') 'REMARK ',tytul(:20),&
1785 ' ENERGY ',etot,' RMS ',rmsd
1789 boxxxshift(1)=int(c(1,nnt)/boxxsize)
1790 ! if (c(1,nnt).lt.0) boxxxshift(1)=boxxxshift(1)-1
1791 boxxxshift(2)=int(c(2,nnt)/boxzsize)
1792 ! if (c(2,nnt).lt.0) boxxxshift(2)=boxxxshift(2)-1
1793 boxxxshift(3)=int(c(3,nnt)/boxzsize)
1794 ! if (c(3,nnt).lt.0) boxxxshift(3)=boxxxshift(3)-1
1803 if (iti.eq.ntyp1_molec(mnum)) then
1806 write (ipdb,'(a)') 'TER'
1813 if ((c(j,i).gt.0).and.(c(j,nnt).lt.0)) then
1814 c(j,i)=dmod(c(j,i),boxxxx(j))+(boxxxshift(j)-1)*boxxxx(j)
1815 else if ((c(j,i).lt.0).and.(c(j,nnt).gt.0)) then
1816 c(j,i)=dmod(c(j,i),boxxxx(j))+(boxxxshift(j)+1)*boxxxx(j)
1818 c(j,i)=dmod(c(j,i),boxxxx(j))+(boxxxshift(j))*boxxxx(j)
1822 write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),&
1823 ires,(c(j,i),j=1,3),1.0d0,tempfac(1,i)
1824 if ((iti.ne.10).and.(mnum.ne.5)) then
1826 write (ipdb,20) iatom,restyp(iti,mnum),chainid(ichain),&
1827 ires,(c(j,nres+i),j=1,3),1.0d0,tempfac(2,i)
1831 write (ipdb,'(a)') 'TER'
1835 if ((itype(i,mnum).eq.ntyp1).or.(mnum.eq.5)) cycle
1836 if (itype(i,mnum).eq.10 .and. itype(i+1,mnum1).ne.ntyp1_molec(mnum1)) then
1837 write (ipdb,30) ica(i),ica(i+1)
1838 else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum1).ne.ntyp1_molec(mnum1)) then
1839 write (ipdb,30) ica(i),ica(i+1),ica(i)+1
1840 else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) then
1841 write (ipdb,30) ica(i),ica(i)+1
1844 if ((itype(nct,molnum(nct)).ne.10).and.(molnum(nct).ne.5)) then
1845 write (ipdb,30) ica(nct),ica(nct)+1
1848 write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
1850 write (ipdb,'(a6)') 'ENDMDL'
1851 10 FORMAT ('ATOM',I7,' CA ',A3,1X,A1,I4,4X,3F8.3,2f6.2)
1852 20 FORMAT ('ATOM',I7,' CB ',A3,1X,A1,I4,4X,3F8.3,2f6.2)
1853 30 FORMAT ('CONECT',8I5)
1855 end subroutine pdboutC
1856 !-----------------------------------------------------------------------------
1857 subroutine cartout(igr,i,etot,free,rmsd,plik)
1858 ! implicit real*8 (a-h,o-z)
1859 ! include 'DIMENSIONS'
1860 ! include 'sizesclu.dat'
1861 ! include 'COMMON.IOUNITS'
1862 ! include 'COMMON.CHAIN'
1863 ! include 'COMMON.VAR'
1864 ! include 'COMMON.LOCAL'
1865 ! include 'COMMON.INTERACT'
1866 ! include 'COMMON.NAMES'
1867 ! include 'COMMON.GEO'
1868 ! include 'COMMON.CLUSTER'
1869 integer :: igr,i,j,k
1870 real(kind=8) :: etot,free,rmsd
1871 character(len=80) :: plik
1872 open (igeom,file=plik,position='append')
1873 write (igeom,'(2e15.5,f10.5,$)') etot,free,rmsd
1874 write (igeom,'(i4,$)') &
1875 nss_all(i),(ihpb_all(j,i),jhpb_all(j,i),j=1,nss_all(i))
1876 write (igeom,'(i10)') iscore(i)
1877 write (igeom,'(8f10.5)') &
1878 ((allcart(k,j,i),k=1,3),j=1,nres),&
1879 ((allcart(k,j+nres,i),k=1,3),j=nnt,nct)
1881 end subroutine cartout
1882 !------------------------------------------------------------------------------
1883 ! subroutine alloc_clust_arrays(n_conf)
1888 ! allocate(diss(maxdist)) !(maxdist)
1889 !el allocate(energy(0:maxconf),totfree(0:maxconf)) !(0:maxconf)
1890 ! allocatable :: enetb !(max_ene,maxstr_proc)
1891 !el allocate(entfac(maxconf)) !(maxconf)
1892 ! allocatable :: totfree_gr !(maxgr)
1893 !el allocate(rcutoff(max_cut+1)) !(max_cut+1)
1895 ! allocatable :: licz,iass !(maxgr)
1896 ! allocatable :: nconf !(maxgr,maxingr)
1897 ! allocatable :: iass_tot !(maxgr,max_cut)
1898 ! allocatable :: list_conf !(maxconf)
1900 !el allocatable :: allcart !(3,maxres2,maxstr_proc)
1901 !el allocate(rmstb(maxconf)) !(maxconf)
1902 !el allocate(mult(nres)) !(maxres)
1903 !el allocatable :: nss_all !(maxstr_proc)
1904 !el allocatable :: ihpb_all,jhpb_all !(maxss,maxstr_proc)
1905 ! allocate(icc(n_conf),iscore(n_conf)) !(maxconf)
1908 ! allocatable :: tempfac !(2,maxres)
1911 !el allocate(beta_h(maxT)) !(maxT)
1913 ! end subroutine alloc_clust_arrays
1914 !-----------------------------------------------------------------------------
1915 !-----------------------------------------------------------------------------