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 io_base, only: rescode
1361 use control, only: setup_var,init_int_table
1362 use conform_compar, only: contact
1364 ! include 'DIMENSIONS'
1365 ! include 'COMMON.IOUNITS'
1366 ! include 'COMMON.GEO'
1367 ! include 'COMMON.VAR'
1368 ! include 'COMMON.INTERACT'
1369 ! include 'COMMON.LOCAL'
1370 ! include 'COMMON.NAMES'
1371 ! include 'COMMON.CHAIN'
1372 ! include 'COMMON.FFIELD'
1373 ! include 'COMMON.SBRIDGE'
1374 ! include 'COMMON.HEADER'
1375 ! include 'COMMON.CONTROL'
1376 ! include 'COMMON.CONTACTS'
1377 ! include 'COMMON.TIME1'
1379 ! include 'COMMON.INFO'
1381 character(len=4) :: sequence(nres) !(maxres)
1382 character(len=800) :: weightcard
1384 real(kind=8) :: x(6*nres) !(maxvar)
1385 integer :: itype_pdb(nres) !(maxres)
1387 integer :: i,j,kkk,mnum
1391 !el allocate(weights(n_ene))
1392 allocate(weights(max_ene))
1393 call alloc_geo_arrays
1394 call alloc_ener_arrays
1395 !-----------------------------
1396 allocate(c(3,2*nres+2)) !(3,maxres2+2) maxres2=2*maxres
1397 allocate(dc(3,0:2*nres+2)) !(3,0:maxres2)
1398 allocate(itype(nres+2,5)) !(maxres)
1399 allocate(molnum(nres+1))
1400 allocate(itel(nres+2))
1414 !--------------------------
1415 ! Read weights of the subsequent energy terms.
1416 call card_concat(weightcard,.true.)
1417 call reada(weightcard,'WSC',wsc,1.0d0)
1418 call reada(weightcard,'WLONG',wsc,wsc)
1419 call reada(weightcard,'WSCP',wscp,1.0d0)
1420 call reada(weightcard,'WELEC',welec,1.0D0)
1421 call reada(weightcard,'WVDWPP',wvdwpp,welec)
1422 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
1423 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
1424 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
1425 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
1426 call reada(weightcard,'WTURN3',wturn3,1.0D0)
1427 call reada(weightcard,'WTURN4',wturn4,1.0D0)
1428 call reada(weightcard,'WTURN6',wturn6,1.0D0)
1429 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
1430 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
1431 call reada(weightcard,'WBOND',wbond,1.0D0)
1432 call reada(weightcard,'WTOR',wtor,1.0D0)
1433 call reada(weightcard,'WTORD',wtor_d,1.0D0)
1434 call reada(weightcard,'WANG',wang,1.0D0)
1435 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
1436 call reada(weightcard,'SCAL14',scal14,0.4D0)
1437 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
1438 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
1439 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
1440 call reada(weightcard,'WCATCAT',wcatcat,0.0d0)
1441 call reada(weightcard,'WCATPROT',wcatprot,0.0d0)
1442 call reada(weightcard,'TEMP0',temp0,300.0d0) !!! el
1443 if (index(weightcard,'SOFT').gt.0) ipot=6
1444 ! 12/1/95 Added weight for the multi-body term WCORR
1445 call reada(weightcard,'WCORRH',wcorr,1.0D0)
1446 if (wcorr4.gt.0.0d0) wcorr=wcorr4
1465 !el weights(19)=wsccor !!!!!!!!!!!!!!!!
1467 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
1468 wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,&
1469 wturn4,wturn6,wsccor
1470 10 format (/'Energy-term weights (unscaled):'// &
1471 'WSCC= ',f10.6,' (SC-SC)'/ &
1472 'WSCP= ',f10.6,' (SC-p)'/ &
1473 'WELEC= ',f10.6,' (p-p electr)'/ &
1474 'WVDWPP= ',f10.6,' (p-p VDW)'/ &
1475 'WBOND= ',f10.6,' (stretching)'/ &
1476 'WANG= ',f10.6,' (bending)'/ &
1477 'WSCLOC= ',f10.6,' (SC local)'/ &
1478 'WTOR= ',f10.6,' (torsional)'/ &
1479 'WTORD= ',f10.6,' (double torsional)'/ &
1480 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ &
1481 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ &
1482 'WCORR4= ',f10.6,' (multi-body 4th order)'/ &
1483 'WCORR5= ',f10.6,' (multi-body 5th order)'/ &
1484 'WCORR6= ',f10.6,' (multi-body 6th order)'/ &
1485 'WTURN3= ',f10.6,' (turns, 3rd order)'/ &
1486 'WTURN4= ',f10.6,' (turns, 4th order)'/ &
1487 'WTURN6= ',f10.6,' (turns, 6th order)'/ &
1488 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
1490 if (wcorr4.gt.0.0d0) then
1491 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',&
1492 'between contact pairs of peptide groups'
1493 write (iout,'(2(a,f5.3/))') &
1494 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,&
1495 'Range of quenching the correlation terms:',2*delt_corr
1496 else if (wcorr.gt.0.0d0) then
1497 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',&
1498 'between contact pairs of peptide groups'
1500 write (iout,'(a,f8.3)') &
1501 'Scaling factor of 1,4 SC-p interactions:',scal14
1502 write (iout,'(a,f8.3)') &
1503 'General scaling factor of SC-p interactions:',scalscp
1504 r0_corr=cutoff_corr-delt_corr
1506 aad(i,1)=scalscp*aad(i,1)
1507 aad(i,2)=scalscp*aad(i,2)
1508 bad(i,1)=scalscp*bad(i,1)
1509 bad(i,2)=scalscp*bad(i,2)
1513 print *,'indpdb=',indpdb,' pdbref=',pdbref
1515 ! Read sequence if not taken from the pdb file.
1516 if (iscode.gt.0) then
1517 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
1519 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
1521 ! Convert sequence to numeric code
1522 do i=1,nres_molec(1)
1525 itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
1527 do i=nres_molec(1)+1,nres_molec(1)+nres_molec(2)
1530 itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
1532 do i=nres_molec(1)+nres_molec(2)+1,nres_molec(1)+nres_molec(2)+nres_molec(5)
1535 itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
1538 print '(20i4)',(itype(i,molnum(i)),i=1,nres)
1542 if (itype(i,molnum(i)).eq.ntyp1_molec(molnum(i))&
1543 .or. itype(i+1,molnum(i+1)).eq.ntyp1_molec(molnum(i+1))) then
1545 if (itype(i,molnum(i)).eq.ntyp1_molec(molnum(i))) then
1549 else if (iabs(itype(i+1,1)).ne.20) then
1551 else if (iabs(itype(i,1)).ne.20) then
1558 write (iout,*) "ITEL"
1560 write (iout,*) i,itype(i,molnum(i)),itel(i)
1563 print *,'Call Read_Bridge.'
1567 print *,'NNT=',NNT,' NCT=',NCT
1568 if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2
1569 if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1
1570 if (nstart.lt.nnt) nstart=nnt
1571 if (nend.gt.nct .or. nend.eq.0) nend=nct
1572 write (iout,*) "nstart",nstart," nend",nend
1575 ! read(inp,'(a)') pdbfile
1576 ! write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
1577 ! open(ipdbin,file=pdbfile,status='old',err=33)
1579 ! 33 write (iout,'(a)') 'Error opening PDB file.'
1582 ! print *,'Begin reading pdb data'
1584 ! print *,'Finished reading pdb data'
1585 ! write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
1587 ! itype_pdb(i)=itype(i)
1590 ! write (iout,'(a,i3)') 'nsup=',nsup
1592 ! if (nsup.le.(nct-nnt+1)) then
1593 ! do i=0,nct-nnt+1-nsup
1594 ! if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1599 ! write (iout,'(a)')
1600 ! & 'Error - sequences to be superposed do not match.'
1603 ! do i=0,nsup-(nct-nnt+1)
1604 ! if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1606 ! nstart_sup=nstart_sup+i
1611 ! write (iout,'(a)')
1612 ! & 'Error - sequences to be superposed do not match.'
1615 ! write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1616 ! & ' nstart_seq=',nstart_seq
1618 write(iout,*)"przed ini_int_tab"
1620 write(iout,*)"po ini_int_tab"
1621 write(iout,*)"przed setup var"
1623 write(iout,*)"po setup var"
1624 write (iout,*) "molread: REFSTR",refstr
1626 if (.not.pdbref) then
1627 call read_angles(inp,*38)
1629 38 write (iout,'(a)') 'Error reading reference structure.'
1631 call mp_stopall(Error_Msg)
1633 stop 'Error reading reference structure'
1642 cref(j,i,kkk)=c(j,i)
1646 call contact(.true.,ncont_ref,icont_ref)
1649 end subroutine molread
1650 !-----------------------------------------------------------------------------
1651 subroutine openunits
1653 ! include 'DIMENSIONS'
1654 use control_data, only: from_cx,from_bx,from_cart
1658 character(len=3) :: liczba
1659 ! include "COMMON.MPI"
1661 ! include 'COMMON.IOUNITS'
1662 ! include 'COMMON.CONTROL'
1663 integer :: lenpre,lenpot !,ilen
1665 character(len=16) :: cformat,cprint
1666 ! character(len=16) ucase
1667 integer :: lenint,lenout
1668 call getenv('INPUT',prefix)
1669 call getenv('OUTPUT',prefout)
1670 call getenv('INTIN',prefintin)
1671 call getenv('COORD',cformat)
1672 call getenv('PRINTCOOR',cprint)
1673 call getenv('SCRATCHDIR',scratchdir)
1676 if (index(ucase(cformat),'CX').gt.0) then
1682 lenout=ilen(prefout)
1683 lenint=ilen(prefintin)
1684 ! Get the names and open the input files
1685 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
1687 write (liczba,'(bz,i3.3)') me
1688 outname=prefout(:lenout)//'_clust.out_'//liczba
1690 outname=prefout(:lenout)//'_clust.out'
1693 intinname=prefintin(:lenint)//'.bx'
1694 else if (from_cx) then
1695 intinname=prefintin(:lenint)//'.cx'
1697 intinname=prefintin(:lenint)//'.int'
1699 rmsname=prefintin(:lenint)//'.rms'
1700 open (jplot,file=prefout(:ilen(prefout))//'.tex',&
1702 open (jrms,file=rmsname,status='unknown')
1703 open(iout,file=outname,status='unknown')
1704 ! Get parameter filenames and open the parameter files.
1705 call getenv('BONDPAR',bondname)
1706 open (ibond,file=bondname,status='old')
1707 call getenv('THETPAR',thetname)
1708 open (ithep,file=thetname,status='old')
1709 call getenv('ROTPAR',rotname)
1710 open (irotam,file=rotname,status='old')
1711 call getenv('TORPAR',torname)
1712 open (itorp,file=torname,status='old')
1713 call getenv('TORDPAR',tordname)
1714 open (itordp,file=tordname,status='old')
1715 call getenv('FOURIER',fouriername)
1716 open (ifourier,file=fouriername,status='old')
1717 call getenv('ELEPAR',elename)
1718 open (ielep,file=elename,status='old')
1719 call getenv('SIDEPAR',sidename)
1720 open (isidep,file=sidename,status='old')
1721 call getenv('SIDEP',sidepname)
1722 open (isidep1,file=sidepname,status="old")
1723 call getenv('SCCORPAR',sccorname)
1724 open (isccor,file=sccorname,status="old")
1725 call getenv('THETPAR_NUCL',thetname_nucl)
1726 open (ithep_nucl,file=thetname_nucl,status='old')
1727 call getenv('ROTPAR_NUCL',rotname_nucl)
1728 open (irotam_nucl,file=rotname_nucl,status='old')
1729 call getenv('TORPAR_NUCL',torname_nucl)
1730 open (itorp_nucl,file=torname_nucl,status='old')
1731 call getenv('TORDPAR_NUCL',tordname_nucl)
1732 open (itordp_nucl,file=tordname_nucl,status='old')
1733 call getenv('SIDEPAR_NUCL',sidename_nucl)
1734 open (isidep_nucl,file=sidename_nucl,status='old')
1735 call getenv('SIDEPAR_SCBASE',sidename_scbase)
1736 open (isidep_scbase,file=sidename_scbase,status='old')
1737 call getenv('PEPPAR_PEPBASE',pepname_pepbase)
1738 open (isidep_pepbase,file=pepname_pepbase,status='old')
1739 call getenv('SCPAR_PHOSPH',pepname_scpho)
1740 open (isidep_scpho,file=pepname_scpho,status='old')
1741 call getenv('PEPPAR_PHOSPH',pepname_peppho)
1742 open (isidep_peppho,file=pepname_peppho,status='old')
1745 call getenv('LIPTRANPAR',liptranname)
1746 open (iliptranpar,file=liptranname,status='old')
1747 call getenv('TUBEPAR',tubename)
1748 open (itube,file=tubename,status='old')
1749 call getenv('IONPAR',ionname)
1750 open (iion,file=ionname,status='old')
1754 ! 8/9/01 In the newest version SCp interaction constants are read from a file
1755 ! Use -DOLDSCP to use hard-coded constants instead.
1757 call getenv('SCPPAR',scpname)
1758 open (iscpp,file=scpname,status='old')
1761 end subroutine openunits
1762 !-----------------------------------------------------------------------------
1764 !-----------------------------------------------------------------------------
1765 subroutine pdboutC(etot,rmsd,tytul)
1767 use energy_data, only: ihpb,jhpb,itype,molnum
1768 ! implicit real*8 (a-h,o-z)
1769 ! include 'DIMENSIONS'
1770 ! include 'COMMON.CONTROL'
1771 ! include 'COMMON.CHAIN'
1772 ! include 'COMMON.INTERACT'
1773 ! include 'COMMON.NAMES'
1774 ! include 'COMMON.IOUNITS'
1775 ! include 'COMMON.HEADER'
1776 ! include 'COMMON.SBRIDGE'
1777 ! include 'COMMON.TEMPFAC'
1778 character(len=50) :: tytul
1779 character(len=1) :: chainid(10)=(/'A','B','C','D','E','F',&
1781 integer :: ica(nres)
1782 real(kind=8) :: etot,rmsd
1783 integer :: iatom,ichain,ires,i,j,iti,mnum,mnum1,boxxxshift(3)
1784 real(kind=8) :: boxxxx(3)
1785 write (ipdb,'(3a,1pe15.5,a,0pf7.2)') 'REMARK ',tytul(:20),&
1786 ' ENERGY ',etot,' RMS ',rmsd
1790 boxxxshift(1)=int(c(1,nnt)/boxxsize)
1791 ! if (c(1,nnt).lt.0) boxxxshift(1)=boxxxshift(1)-1
1792 boxxxshift(2)=int(c(2,nnt)/boxzsize)
1793 ! if (c(2,nnt).lt.0) boxxxshift(2)=boxxxshift(2)-1
1794 boxxxshift(3)=int(c(3,nnt)/boxzsize)
1795 ! if (c(3,nnt).lt.0) boxxxshift(3)=boxxxshift(3)-1
1804 if (iti.eq.ntyp1_molec(mnum)) then
1807 write (ipdb,'(a)') 'TER'
1814 if ((c(j,i).gt.0).and.(c(j,nnt).lt.0)) then
1815 c(j,i)=dmod(c(j,i),boxxxx(j))+(boxxxshift(j)-1)*boxxxx(j)
1816 else if ((c(j,i).lt.0).and.(c(j,nnt).gt.0)) then
1817 c(j,i)=dmod(c(j,i),boxxxx(j))+(boxxxshift(j)+1)*boxxxx(j)
1819 c(j,i)=dmod(c(j,i),boxxxx(j))+(boxxxshift(j))*boxxxx(j)
1823 write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),&
1824 ires,(c(j,i),j=1,3),1.0d0,tempfac(1,i)
1825 if ((iti.ne.10).and.(mnum.ne.5)) then
1827 write (ipdb,20) iatom,restyp(iti,mnum),chainid(ichain),&
1828 ires,(c(j,nres+i),j=1,3),1.0d0,tempfac(2,i)
1832 write (ipdb,'(a)') 'TER'
1836 if ((itype(i,mnum).eq.ntyp1).or.(mnum.eq.5)) cycle
1837 if (itype(i,mnum).eq.10 .and. itype(i+1,mnum1).ne.ntyp1_molec(mnum1)) then
1838 write (ipdb,30) ica(i),ica(i+1)
1839 else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum1).ne.ntyp1_molec(mnum1)) then
1840 write (ipdb,30) ica(i),ica(i+1),ica(i)+1
1841 else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) then
1842 write (ipdb,30) ica(i),ica(i)+1
1845 if ((itype(nct,molnum(nct)).ne.10).and.(molnum(nct).ne.5)) then
1846 write (ipdb,30) ica(nct),ica(nct)+1
1849 write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
1851 write (ipdb,'(a6)') 'ENDMDL'
1852 10 FORMAT ('ATOM',I7,' CA ',A3,1X,A1,I4,4X,3F8.3,2f6.2)
1853 20 FORMAT ('ATOM',I7,' CB ',A3,1X,A1,I4,4X,3F8.3,2f6.2)
1854 30 FORMAT ('CONECT',8I5)
1856 end subroutine pdboutC
1857 !-----------------------------------------------------------------------------
1858 subroutine cartout(igr,i,etot,free,rmsd,plik)
1859 ! implicit real*8 (a-h,o-z)
1860 ! include 'DIMENSIONS'
1861 ! include 'sizesclu.dat'
1862 ! include 'COMMON.IOUNITS'
1863 ! include 'COMMON.CHAIN'
1864 ! include 'COMMON.VAR'
1865 ! include 'COMMON.LOCAL'
1866 ! include 'COMMON.INTERACT'
1867 ! include 'COMMON.NAMES'
1868 ! include 'COMMON.GEO'
1869 ! include 'COMMON.CLUSTER'
1870 integer :: igr,i,j,k
1871 real(kind=8) :: etot,free,rmsd
1872 character(len=80) :: plik
1873 open (igeom,file=plik,position='append')
1874 write (igeom,'(2e15.5,f10.5,$)') etot,free,rmsd
1875 write (igeom,'(i4,$)') &
1876 nss_all(i),(ihpb_all(j,i),jhpb_all(j,i),j=1,nss_all(i))
1877 write (igeom,'(i10)') iscore(i)
1878 write (igeom,'(8f10.5)') &
1879 ((allcart(k,j,i),k=1,3),j=1,nres),&
1880 ((allcart(k,j+nres,i),k=1,3),j=nnt,nct)
1882 end subroutine cartout
1883 !------------------------------------------------------------------------------
1884 ! subroutine alloc_clust_arrays(n_conf)
1889 ! allocate(diss(maxdist)) !(maxdist)
1890 !el allocate(energy(0:maxconf),totfree(0:maxconf)) !(0:maxconf)
1891 ! allocatable :: enetb !(max_ene,maxstr_proc)
1892 !el allocate(entfac(maxconf)) !(maxconf)
1893 ! allocatable :: totfree_gr !(maxgr)
1894 !el allocate(rcutoff(max_cut+1)) !(max_cut+1)
1896 ! allocatable :: licz,iass !(maxgr)
1897 ! allocatable :: nconf !(maxgr,maxingr)
1898 ! allocatable :: iass_tot !(maxgr,max_cut)
1899 ! allocatable :: list_conf !(maxconf)
1901 !el allocatable :: allcart !(3,maxres2,maxstr_proc)
1902 !el allocate(rmstb(maxconf)) !(maxconf)
1903 !el allocate(mult(nres)) !(maxres)
1904 !el allocatable :: nss_all !(maxstr_proc)
1905 !el allocatable :: ihpb_all,jhpb_all !(maxss,maxstr_proc)
1906 ! allocate(icc(n_conf),iscore(n_conf)) !(maxconf)
1909 ! allocatable :: tempfac !(2,maxres)
1912 !el allocate(beta_h(maxT)) !(maxT)
1914 ! end subroutine alloc_clust_arrays
1915 !-----------------------------------------------------------------------------
1916 !-----------------------------------------------------------------------------