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.)
892 write (iout,*) "Check atom",j
894 if (itype(j,mnum).eq.ntyp1_molec(mnum)) cycle
895 if (itype(j+1,molnum(j+1)).eq.ntyp1_molec(molnum(j+1))) cycle
897 if ((vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0).and.(mnum.ne.5)) then
899 if (itel(j).ne.0 .and. itel(j-1).ne.0) then
900 write (iout,*) "Conformation",jjj,jj+1
901 write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),itel(j),&
903 write (iout,*) "The Cartesian geometry is:"
904 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
905 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
906 write (iout,*) "The internal geometry is:"
907 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
908 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
909 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
910 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
911 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
912 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
914 "This conformation WILL NOT be added to the database."
924 if (itype(j,1).ne.10 .and. (itype(j,mnum).ne.ntyp1_molec(mnum)) &
925 .and. (vbld(nres+j)-dsc(iabs(itj))) &
927 write (iout,*) "Conformation",jjj,jj+1
928 write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j)
929 write (iout,*) "The Cartesian geometry is:"
930 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
931 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
932 write (iout,*) "The internal geometry is:"
933 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
934 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
935 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
936 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
937 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
938 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
940 "This conformation WILL NOT be added to the database."
945 if (theta(j).le.0.0d0) then
947 "Zero theta angle(s) in conformation",jjj,jj+1
948 write (iout,*) "The Cartesian geometry is:"
949 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
950 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
951 write (iout,*) "The internal geometry is:"
952 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
953 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
954 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
955 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
956 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
957 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
959 "This conformation WILL NOT be added to the database."
962 if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad
966 write (iout,*) "Conformation",jjj,jj
967 write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
968 write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
969 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
970 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
971 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
972 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
973 write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct)
974 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
975 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
976 write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
977 write (iout,'(e15.5,16i5)') entfac(icount+1)
978 ! & iscore(icount+1,0)
981 call store_cconf_from_file(jj,icount)
982 if (icount.eq.maxstr_proc) then
984 write (iout,* ) "jj_old",jj_old," jj",jj
986 call write_and_send_cconf(icount,jj_old,jj,Next)
991 end subroutine add_new_cconf
992 !------------------------------------------------------------------------------
993 subroutine store_cconf_from_file(jj,icount)
995 use energy_data, only: ihpb,jhpb
997 ! include "DIMENSIONS"
998 ! include "sizesclu.dat"
999 ! include "COMMON.CLUSTER"
1000 ! include "COMMON.CHAIN"
1001 ! include "COMMON.SBRIDGE"
1002 ! include "COMMON.INTERACT"
1003 ! include "COMMON.IOUNITS"
1004 ! include "COMMON.VAR"
1005 integer :: i,j,jj,icount
1006 ! Store the conformation that has been read in
1009 allcart(j,i,icount)=c(j,i)
1014 ihpb_all(i,icount)=ihpb(i)
1015 jhpb_all(i,icount)=jhpb(i)
1018 end subroutine store_cconf_from_file
1019 !------------------------------------------------------------------------------
1020 subroutine write_and_send_cconf(icount,jj_old,jj,Next)
1023 ! include "DIMENSIONS"
1024 ! include "sizesclu.dat"
1029 ! include "COMMON.MPI"
1031 ! include "COMMON.CHAIN"
1032 ! include "COMMON.SBRIDGE"
1033 ! include "COMMON.INTERACT"
1034 ! include "COMMON.IOUNITS"
1035 ! include "COMMON.CLUSTER"
1036 ! include "COMMON.VAR"
1037 integer :: icount,jj_old,jj,Next
1038 ! Write the structures to a scratch file
1040 ! Master sends the portion of conformations that have been read in to the neighbor
1042 write (iout,*) "Processor",me," entered WRITE_AND_SEND_CONF"
1045 call MPI_Send(icount,1,MPI_INTEGER,Next,570,MPI_COMM_WORLD,IERROR)
1046 call MPI_Send(nss_all(1),icount,MPI_INTEGER,&
1047 Next,571,MPI_COMM_WORLD,IERROR)
1048 call MPI_Send(ihpb_all(1,1),icount,MPI_INTEGER,&
1049 Next,572,MPI_COMM_WORLD,IERROR)
1050 call MPI_Send(jhpb_all(1,1),icount,MPI_INTEGER,&
1051 Next,573,MPI_COMM_WORLD,IERROR)
1052 call MPI_Send(rmstb(jj_old),icount,MPI_DOUBLE_PRECISION,&
1053 Next,577,MPI_COMM_WORLD,IERROR)
1054 call MPI_Send(entfac(jj_old),icount,MPI_DOUBLE_PRECISION,&
1055 Next,579,MPI_COMM_WORLD,IERROR)
1056 call MPI_Send(allcart(1,1,1),3*icount*2*nres,&
1057 MPI_REAL,Next,580,MPI_COMM_WORLD,IERROR)
1059 call dawrite_ccoords(jj_old,jj,icbase)
1061 end subroutine write_and_send_cconf
1062 !------------------------------------------------------------------------------
1064 subroutine receive_and_pass_cconf(icount,jj_old,jj,Previous,Next)
1068 ! include "DIMENSIONS"
1069 ! include "sizesclu.dat"
1071 integer :: IERROR !,STATUS(MPI_STATUS_SIZE)
1072 ! include "COMMON.MPI"
1073 ! include "COMMON.CHAIN"
1074 ! include "COMMON.SBRIDGE"
1075 ! include "COMMON.INTERACT"
1076 ! include "COMMON.IOUNITS"
1077 ! include "COMMON.VAR"
1078 ! include "COMMON.GEO"
1079 ! include "COMMON.CLUSTER"
1080 integer :: i,j,k,icount,jj_old,jj,Previous,Next
1083 write (iout,*) "Processor",me," entered RECEIVE_AND_PASS_CONF"
1086 do while (icount.gt.0)
1087 call MPI_Recv(icount,1,MPI_INTEGER,Previous,570,MPI_COMM_WORLD,&
1089 call MPI_Send(icount,1,MPI_INTEGER,Next,570,MPI_COMM_WORLD,&
1092 write (iout,*) "Processor",me," icount",icount
1094 if (icount.eq.0) return
1095 call MPI_Recv(nss_all(1),icount,MPI_INTEGER,&
1096 Previous,571,MPI_COMM_WORLD,STATUS,IERROR)
1097 call MPI_Send(nss_all(1),icount,MPI_INTEGER,&
1098 Next,571,MPI_COMM_WORLD,IERROR)
1099 call MPI_Recv(ihpb_all(1,1),icount,MPI_INTEGER,&
1100 Previous,572,MPI_COMM_WORLD,STATUS,IERROR)
1101 call MPI_Send(ihpb_all(1,1),icount,MPI_INTEGER,&
1102 Next,572,MPI_COMM_WORLD,IERROR)
1103 call MPI_Recv(jhpb_all(1,1),icount,MPI_INTEGER,&
1104 Previous,573,MPI_COMM_WORLD,STATUS,IERROR)
1105 call MPI_Send(jhpb_all(1,1),icount,MPI_INTEGER,&
1106 Next,573,MPI_COMM_WORLD,IERROR)
1107 call MPI_Recv(rmstb(jj_old),icount,MPI_DOUBLE_PRECISION,&
1108 Previous,577,MPI_COMM_WORLD,STATUS,IERROR)
1109 call MPI_Send(rmstb(jj_old),icount,MPI_DOUBLE_PRECISION,&
1110 Next,577,MPI_COMM_WORLD,IERROR)
1111 call MPI_Recv(entfac(jj_old),icount,MPI_DOUBLE_PRECISION,&
1112 Previous,579,MPI_COMM_WORLD,STATUS,IERROR)
1113 call MPI_Send(entfac(jj_old),icount,MPI_DOUBLE_PRECISION,&
1114 Next,579,MPI_COMM_WORLD,IERROR)
1115 call MPI_Recv(allcart(1,1,1),3*icount*2*nres,&
1116 MPI_REAL,Previous,580,MPI_COMM_WORLD,STATUS,IERROR)
1117 call MPI_Send(allcart(1,1,1),3*icount*2*nres,&
1118 MPI_REAL,Next,580,MPI_COMM_WORLD,IERROR)
1120 call dawrite_ccoords(jj_old,jj,icbase)
1123 write (iout,*) "Processor",me," received",icount," conformations"
1125 write (iout,'(8f10.4)') (allcart(l,k,i),l=1,3,k=1,nres)
1126 write (iout,'(8f10.4)')((allcart(l,k,i+nres),l=1,3,k=nnt,nct)
1127 write (iout,'(e15.5,16i5)') entfac(i)
1132 end subroutine receive_and_pass_cconf
1134 !------------------------------------------------------------------------------
1135 subroutine daread_ccoords(istart_conf,iend_conf)
1138 ! include "DIMENSIONS"
1139 ! include "sizesclu.dat"
1143 ! include "COMMON.MPI"
1145 ! include "COMMON.CHAIN"
1146 ! include "COMMON.CLUSTER"
1147 ! include "COMMON.IOUNITS"
1148 ! include "COMMON.INTERACT"
1149 ! include "COMMON.VAR"
1150 ! include "COMMON.SBRIDGE"
1151 ! include "COMMON.GEO"
1152 integer :: istart_conf,iend_conf
1153 integer :: i,j,ij,ii,iii
1155 character(len=16) :: form,acc
1156 character(len=32) :: nam
1158 ! Read conformations off a DA scratchfile.
1161 write (iout,*) "DAREAD_COORDS"
1162 write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
1163 inquire(unit=icbase,name=nam,recl=len,form=form,access=acc)
1164 write (iout,*) "len=",len," form=",form," acc=",acc
1165 write (iout,*) "nam=",nam
1168 do ii=istart_conf,iend_conf
1169 ij = ii - istart_conf + 1
1172 write (iout,*) "Reading binary file, record",iii," ii",ii
1175 read(icbase,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres),&
1176 ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
1177 nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss),&
1178 entfac(ii),rmstb(ii)
1180 write (iout,*) ii,iii,ij,entfac(ii)
1181 write (iout,'(8f10.5)') ((allcart(j,i,ij),j=1,3),i=1,nres)
1182 write (iout,'(8f10.4)') ((allcart(j,i,ij),j=1,3),&
1183 i=nnt+nres,nct+nres)
1184 write (iout,'(2e15.5)') entfac(ij)
1185 write (iout,'(16i5)') nss_all(ij),(ihpb_all(i,ij),&
1186 jhpb_all(i,ij),i=1,nss)
1191 end subroutine daread_ccoords
1192 !------------------------------------------------------------------------------
1193 subroutine dawrite_ccoords(istart_conf,iend_conf,unit_out)
1196 ! include "DIMENSIONS"
1197 ! include "sizesclu.dat"
1201 ! include "COMMON.MPI"
1203 ! include "COMMON.CHAIN"
1204 ! include "COMMON.INTERACT"
1205 ! include "COMMON.IOUNITS"
1206 ! include "COMMON.VAR"
1207 ! include "COMMON.SBRIDGE"
1208 ! include "COMMON.GEO"
1209 ! include "COMMON.CLUSTER"
1210 integer :: istart_conf,iend_conf
1211 integer :: i,j,ii,ij,iii,unit_out
1213 character(len=16) :: form,acc
1214 character(len=32) :: nam
1216 ! Write conformations to a DA scratchfile.
1219 write (iout,*) "DAWRITE_COORDS"
1220 write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
1221 write (iout,*) "lenrec",lenrec
1222 inquire(unit=unit_out,name=nam,recl=len,form=form,access=acc)
1223 write (iout,*) "len=",len," form=",form," acc=",acc
1224 write (iout,*) "nam=",nam
1227 do ii=istart_conf,iend_conf
1229 ij = ii - istart_conf + 1
1231 write (iout,*) "Writing binary file, record",iii," ii",ii
1234 write(unit_out,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres),&
1235 ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
1236 nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss_all(ij)),&
1237 entfac(ii),rmstb(ii)
1239 write (iout,'(8f10.5)') ((allcart(j,i,ij),j=1,3),i=1,nres)
1240 write (iout,'(8f10.4)') ((allcart(j,i,ij),j=1,3),i=nnt+nres,&
1242 write (iout,'(2e15.5)') entfac(ij)
1243 write (iout,'(16i5)') nss_all(ij),(ihpb(i,ij),jhpb(i,ij),i=1,&
1249 end subroutine dawrite_ccoords
1250 !-----------------------------------------------------------------------------
1252 !-----------------------------------------------------------------------------
1253 subroutine read_control
1255 ! Read molecular data
1257 use energy_data, only: rescale_mode,distchainmax,ipot !,temp0
1258 use control_data, only: titel,outpdb,outmol2,refstr,pdbref,&
1259 iscode,symetr,punch_dist,print_dist,nstart,nend,&
1260 caonly,iopt,efree,lprint_cart,lprint_int,rlamb_ele,&
1261 r_cut_ele,nclust,tor_mode,scelemode
1263 ! include 'DIMENSIONS'
1264 ! include 'sizesclu.dat'
1265 ! include 'COMMON.IOUNITS'
1266 ! include 'COMMON.TIME1'
1267 ! include 'COMMON.SBRIDGE'
1268 ! include 'COMMON.CONTROL'
1269 ! include 'COMMON.CLUSTER'
1270 ! include 'COMMON.CHAIN'
1271 ! include 'COMMON.HEADER'
1272 ! include 'COMMON.FFIELD'
1273 ! include 'COMMON.FREE'
1274 ! include 'COMMON.INTERACT'
1275 character(len=320) :: controlcard !,ucase
1277 ! include 'COMMON.INFO'
1281 read (INP,'(a80)') titel
1282 call card_concat(controlcard,.true.)
1284 call readi(controlcard,'NRES',nres_molec(1),0)
1285 call readi(controlcard,"NRES_NUCL",nres_molec(2),0)
1286 call readi(controlcard,"NRES_CAT",nres_molec(5),0)
1289 nres=nres_molec(i)+nres
1292 ! call alloc_clust_arrays
1293 allocate(rcutoff(max_cut+1)) !(max_cut+1)
1294 allocate(beta_h(maxT)) !(maxT)
1295 allocate(mult(nres)) !(maxres)
1298 call readi(controlcard,'RESCALE',rescale_mode,2)
1299 call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
1300 write (iout,*) "DISTCHAINMAX",distchainmax
1301 call readi(controlcard,'PDBOUT',outpdb,0)
1302 call readi(controlcard,'MOL2OUT',outmol2,0)
1303 refstr=(index(controlcard,'REFSTR').gt.0)
1304 write (iout,*) "REFSTR",refstr
1305 pdbref=(index(controlcard,'PDBREF').gt.0)
1306 iscode=index(controlcard,'ONE_LETTER')
1307 tree=(index(controlcard,'MAKE_TREE').gt.0)
1308 min_var=(index(controlcard,'MINVAR').gt.0)
1309 plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
1310 punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
1311 call readi(controlcard,'NCUT',ncut,1)
1312 call readi(controlcard,'SYM',symetr,1)
1313 write (iout,*) 'sym', symetr
1314 call readi(controlcard,'NSTART',nstart,0)
1315 call reada(controlcard,'BOXX',boxxsize,100.0d0)
1316 call reada(controlcard,'BOXY',boxysize,100.0d0)
1317 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
1318 call readi(controlcard,'NCLUST',nclust,5)
1319 ! ions=index(controlcard,"IONS").gt.0
1320 call readi(controlcard,"SCELEMODE",scelemode,0)
1321 print *,"SCELE",scelemode
1322 call readi(controlcard,'TORMODE',tor_mode,0)
1323 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1324 write(iout,*) "torsional and valence angle mode",tor_mode
1325 call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
1326 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
1327 write(iout,*) "R_CUT_ELE=",r_cut_ele
1328 call readi(controlcard,'NEND',nend,0)
1329 call reada(controlcard,'ECUT',ecut,10.0d0)
1330 call reada(controlcard,'PROB',prob_limit,0.99d0)
1331 write (iout,*) "Probability limit",prob_limit
1332 lgrp=(index(controlcard,'LGRP').gt.0)
1333 caonly=(index(controlcard,'CA_ONLY').gt.0)
1334 print_dist=(index(controlcard,'PRINT_DIST').gt.0)
1335 call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
1336 call readi(controlcard,'IOPT',iopt,2)
1337 lside = index(controlcard,"SIDE").gt.0
1338 efree = index(controlcard,"EFREE").gt.0
1339 call readi(controlcard,'NTEMP',nT,1)
1340 write (iout,*) "nT",nT
1341 !el call reada(controlcard,'TEMP0',temp0,300.0d0) !el
1342 call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
1343 write (iout,*) "nT",nT
1344 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
1346 beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
1348 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
1349 lprint_cart=index(controlcard,"PRINT_CART") .gt.0
1350 lprint_int=index(controlcard,"PRINT_INT") .gt.0
1353 end subroutine read_control
1354 !-----------------------------------------------------------------------------
1357 ! Read molecular data.
1359 use geometry_data, only: nsup,cref,nres0,nstart_sup,nstart_seq,dc
1360 use energy_data!, only: wsc,wscp,welec,wbond,wstrain,wtor,wtor_d,&
1361 ! wang,wscloc,wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,&
1362 ! wturn3,wturn4,wturn6,wvdwpp,weights
1363 use control_data, only: titel,nstart,nend,pdbref,refstr,iscode,&
1365 use geometry, only: chainbuild,alloc_geo_arrays
1366 use energy, only: alloc_ener_arrays
1367 use io_base, only: rescode
1368 use control, only: setup_var,init_int_table
1369 use conform_compar, only: contact
1371 ! include 'DIMENSIONS'
1372 ! include 'COMMON.IOUNITS'
1373 ! include 'COMMON.GEO'
1374 ! include 'COMMON.VAR'
1375 ! include 'COMMON.INTERACT'
1376 ! include 'COMMON.LOCAL'
1377 ! include 'COMMON.NAMES'
1378 ! include 'COMMON.CHAIN'
1379 ! include 'COMMON.FFIELD'
1380 ! include 'COMMON.SBRIDGE'
1381 ! include 'COMMON.HEADER'
1382 ! include 'COMMON.CONTROL'
1383 ! include 'COMMON.CONTACTS'
1384 ! include 'COMMON.TIME1'
1386 ! include 'COMMON.INFO'
1388 character(len=4) :: sequence(nres) !(maxres)
1389 character(len=800) :: weightcard
1391 real(kind=8) :: x(6*nres) !(maxvar)
1392 integer :: itype_pdb(nres) !(maxres)
1394 integer :: i,j,kkk,mnum
1398 !el allocate(weights(n_ene))
1399 allocate(weights(max_ene))
1400 call alloc_geo_arrays
1401 call alloc_ener_arrays
1402 !-----------------------------
1403 allocate(c(3,2*nres+2)) !(3,maxres2+2) maxres2=2*maxres
1404 allocate(dc(3,0:2*nres+2)) !(3,0:maxres2)
1405 allocate(itype(nres+2,5)) !(maxres)
1406 allocate(molnum(nres+1))
1407 allocate(itel(nres+2))
1421 !--------------------------
1422 ! Read weights of the subsequent energy terms.
1423 call card_concat(weightcard,.true.)
1424 call reada(weightcard,'WSC',wsc,1.0d0)
1425 call reada(weightcard,'WLONG',wsc,wsc)
1426 call reada(weightcard,'WSCP',wscp,1.0d0)
1427 call reada(weightcard,'WELEC',welec,1.0D0)
1428 call reada(weightcard,'WVDWPP',wvdwpp,welec)
1429 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
1430 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
1431 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
1432 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
1433 call reada(weightcard,'WTURN3',wturn3,1.0D0)
1434 call reada(weightcard,'WTURN4',wturn4,1.0D0)
1435 call reada(weightcard,'WTURN6',wturn6,1.0D0)
1436 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
1437 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
1438 call reada(weightcard,'WBOND',wbond,1.0D0)
1439 call reada(weightcard,'WTOR',wtor,1.0D0)
1440 call reada(weightcard,'WTORD',wtor_d,1.0D0)
1441 call reada(weightcard,'WANG',wang,1.0D0)
1442 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
1443 call reada(weightcard,'SCAL14',scal14,0.4D0)
1444 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
1445 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
1446 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
1447 call reada(weightcard,'WCATCAT',wcatcat,0.0d0)
1448 call reada(weightcard,'WCATPROT',wcatprot,0.0d0)
1449 call reada(weightcard,'TEMP0',temp0,300.0d0) !!! el
1450 if (index(weightcard,'SOFT').gt.0) ipot=6
1451 ! 12/1/95 Added weight for the multi-body term WCORR
1452 call reada(weightcard,'WCORRH',wcorr,1.0D0)
1453 if (wcorr4.gt.0.0d0) wcorr=wcorr4
1472 !el weights(19)=wsccor !!!!!!!!!!!!!!!!
1474 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
1475 wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,&
1476 wturn4,wturn6,wsccor
1477 10 format (/'Energy-term weights (unscaled):'// &
1478 'WSCC= ',f10.6,' (SC-SC)'/ &
1479 'WSCP= ',f10.6,' (SC-p)'/ &
1480 'WELEC= ',f10.6,' (p-p electr)'/ &
1481 'WVDWPP= ',f10.6,' (p-p VDW)'/ &
1482 'WBOND= ',f10.6,' (stretching)'/ &
1483 'WANG= ',f10.6,' (bending)'/ &
1484 'WSCLOC= ',f10.6,' (SC local)'/ &
1485 'WTOR= ',f10.6,' (torsional)'/ &
1486 'WTORD= ',f10.6,' (double torsional)'/ &
1487 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ &
1488 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ &
1489 'WCORR4= ',f10.6,' (multi-body 4th order)'/ &
1490 'WCORR5= ',f10.6,' (multi-body 5th order)'/ &
1491 'WCORR6= ',f10.6,' (multi-body 6th order)'/ &
1492 'WTURN3= ',f10.6,' (turns, 3rd order)'/ &
1493 'WTURN4= ',f10.6,' (turns, 4th order)'/ &
1494 'WTURN6= ',f10.6,' (turns, 6th order)'/ &
1495 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
1497 if (wcorr4.gt.0.0d0) then
1498 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',&
1499 'between contact pairs of peptide groups'
1500 write (iout,'(2(a,f5.3/))') &
1501 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,&
1502 'Range of quenching the correlation terms:',2*delt_corr
1503 else if (wcorr.gt.0.0d0) then
1504 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',&
1505 'between contact pairs of peptide groups'
1507 write (iout,'(a,f8.3)') &
1508 'Scaling factor of 1,4 SC-p interactions:',scal14
1509 write (iout,'(a,f8.3)') &
1510 'General scaling factor of SC-p interactions:',scalscp
1511 r0_corr=cutoff_corr-delt_corr
1513 aad(i,1)=scalscp*aad(i,1)
1514 aad(i,2)=scalscp*aad(i,2)
1515 bad(i,1)=scalscp*bad(i,1)
1516 bad(i,2)=scalscp*bad(i,2)
1520 print *,'indpdb=',indpdb,' pdbref=',pdbref
1522 ! Read sequence if not taken from the pdb file.
1523 if (iscode.gt.0) then
1524 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
1526 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
1528 ! Convert sequence to numeric code
1529 do i=1,nres_molec(1)
1532 itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
1534 do i=nres_molec(1)+1,nres_molec(1)+nres_molec(2)
1537 itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
1539 do i=nres_molec(1)+nres_molec(2)+1,nres_molec(1)+nres_molec(2)+nres_molec(5)
1542 itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
1545 print '(20i4)',(itype(i,molnum(i)),i=1,nres)
1549 if (itype(i,molnum(i)).eq.ntyp1_molec(molnum(i))&
1550 .or. itype(i+1,molnum(i+1)).eq.ntyp1_molec(molnum(i+1))) then
1552 if (itype(i,molnum(i)).eq.ntyp1_molec(molnum(i))) then
1556 else if (iabs(itype(i+1,1)).ne.20) then
1558 else if (iabs(itype(i,1)).ne.20) then
1565 write (iout,*) "ITEL"
1567 write (iout,*) i,itype(i,molnum(i)),itel(i)
1570 print *,'Call Read_Bridge.'
1574 print *,'NNT=',NNT,' NCT=',NCT
1575 if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2
1576 if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1
1577 if (nstart.lt.nnt) nstart=nnt
1578 if (nend.gt.nct .or. nend.eq.0) nend=nct
1579 write (iout,*) "nstart",nstart," nend",nend
1582 ! read(inp,'(a)') pdbfile
1583 ! write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
1584 ! open(ipdbin,file=pdbfile,status='old',err=33)
1586 ! 33 write (iout,'(a)') 'Error opening PDB file.'
1589 ! print *,'Begin reading pdb data'
1591 ! print *,'Finished reading pdb data'
1592 ! write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
1594 ! itype_pdb(i)=itype(i)
1597 ! write (iout,'(a,i3)') 'nsup=',nsup
1599 ! if (nsup.le.(nct-nnt+1)) then
1600 ! do i=0,nct-nnt+1-nsup
1601 ! if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1606 ! write (iout,'(a)')
1607 ! & 'Error - sequences to be superposed do not match.'
1610 ! do i=0,nsup-(nct-nnt+1)
1611 ! if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1613 ! nstart_sup=nstart_sup+i
1618 ! write (iout,'(a)')
1619 ! & 'Error - sequences to be superposed do not match.'
1622 ! write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1623 ! & ' nstart_seq=',nstart_seq
1625 write(iout,*)"przed ini_int_tab"
1627 write(iout,*)"po ini_int_tab"
1628 write(iout,*)"przed setup var"
1630 write(iout,*)"po setup var"
1631 write (iout,*) "molread: REFSTR",refstr
1633 if (.not.pdbref) then
1634 call read_angles(inp,*38)
1636 38 write (iout,'(a)') 'Error reading reference structure.'
1638 call mp_stopall(Error_Msg)
1640 stop 'Error reading reference structure'
1649 cref(j,i,kkk)=c(j,i)
1653 call contact(.true.,ncont_ref,icont_ref)
1656 end subroutine molread
1657 !-----------------------------------------------------------------------------
1658 subroutine openunits
1660 ! include 'DIMENSIONS'
1661 use control_data, only: from_cx,from_bx,from_cart
1665 character(len=3) :: liczba
1666 ! include "COMMON.MPI"
1668 ! include 'COMMON.IOUNITS'
1669 ! include 'COMMON.CONTROL'
1670 integer :: lenpre,lenpot !,ilen
1672 character(len=16) :: cformat,cprint
1673 ! character(len=16) ucase
1674 integer :: lenint,lenout
1675 call getenv('INPUT',prefix)
1676 call getenv('OUTPUT',prefout)
1677 call getenv('INTIN',prefintin)
1678 call getenv('COORD',cformat)
1679 call getenv('PRINTCOOR',cprint)
1680 call getenv('SCRATCHDIR',scratchdir)
1683 if (index(ucase(cformat),'CX').gt.0) then
1689 lenout=ilen(prefout)
1690 lenint=ilen(prefintin)
1691 ! Get the names and open the input files
1692 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
1694 write (liczba,'(bz,i3.3)') me
1695 outname=prefout(:lenout)//'_clust.out_'//liczba
1697 outname=prefout(:lenout)//'_clust.out'
1700 intinname=prefintin(:lenint)//'.bx'
1701 else if (from_cx) then
1702 intinname=prefintin(:lenint)//'.cx'
1704 intinname=prefintin(:lenint)//'.int'
1706 rmsname=prefintin(:lenint)//'.rms'
1707 open (jplot,file=prefout(:ilen(prefout))//'.tex',&
1709 open (jrms,file=rmsname,status='unknown')
1710 open(iout,file=outname,status='unknown')
1711 ! Get parameter filenames and open the parameter files.
1712 call getenv('BONDPAR',bondname)
1713 open (ibond,file=bondname,status='old')
1714 call getenv('THETPAR',thetname)
1715 open (ithep,file=thetname,status='old')
1716 call getenv('ROTPAR',rotname)
1717 open (irotam,file=rotname,status='old')
1718 call getenv('TORPAR',torname)
1719 open (itorp,file=torname,status='old')
1720 call getenv('TORDPAR',tordname)
1721 open (itordp,file=tordname,status='old')
1722 call getenv('FOURIER',fouriername)
1723 open (ifourier,file=fouriername,status='old')
1724 call getenv('ELEPAR',elename)
1725 open (ielep,file=elename,status='old')
1726 call getenv('SIDEPAR',sidename)
1727 open (isidep,file=sidename,status='old')
1728 call getenv('SIDEP',sidepname)
1729 open (isidep1,file=sidepname,status="old")
1730 call getenv('SCCORPAR',sccorname)
1731 open (isccor,file=sccorname,status="old")
1732 call getenv('THETPAR_NUCL',thetname_nucl)
1733 open (ithep_nucl,file=thetname_nucl,status='old')
1734 call getenv('ROTPAR_NUCL',rotname_nucl)
1735 open (irotam_nucl,file=rotname_nucl,status='old')
1736 call getenv('TORPAR_NUCL',torname_nucl)
1737 open (itorp_nucl,file=torname_nucl,status='old')
1738 call getenv('TORDPAR_NUCL',tordname_nucl)
1739 open (itordp_nucl,file=tordname_nucl,status='old')
1740 call getenv('SIDEPAR_NUCL',sidename_nucl)
1741 open (isidep_nucl,file=sidename_nucl,status='old')
1742 call getenv('SIDEPAR_SCBASE',sidename_scbase)
1743 open (isidep_scbase,file=sidename_scbase,status='old')
1744 call getenv('PEPPAR_PEPBASE',pepname_pepbase)
1745 open (isidep_pepbase,file=pepname_pepbase,status='old')
1746 call getenv('SCPAR_PHOSPH',pepname_scpho)
1747 open (isidep_scpho,file=pepname_scpho,status='old')
1748 call getenv('PEPPAR_PHOSPH',pepname_peppho)
1749 open (isidep_peppho,file=pepname_peppho,status='old')
1752 call getenv('LIPTRANPAR',liptranname)
1753 open (iliptranpar,file=liptranname,status='old')
1754 call getenv('TUBEPAR',tubename)
1755 open (itube,file=tubename,status='old')
1756 call getenv('IONPAR',ionname)
1757 open (iion,file=ionname,status='old')
1761 ! 8/9/01 In the newest version SCp interaction constants are read from a file
1762 ! Use -DOLDSCP to use hard-coded constants instead.
1764 call getenv('SCPPAR',scpname)
1765 open (iscpp,file=scpname,status='old')
1768 end subroutine openunits
1769 !-----------------------------------------------------------------------------
1771 !-----------------------------------------------------------------------------
1772 subroutine pdboutC(etot,rmsd,tytul)
1774 use energy_data, only: ihpb,jhpb,itype,molnum
1775 ! implicit real*8 (a-h,o-z)
1776 ! include 'DIMENSIONS'
1777 ! include 'COMMON.CONTROL'
1778 ! include 'COMMON.CHAIN'
1779 ! include 'COMMON.INTERACT'
1780 ! include 'COMMON.NAMES'
1781 ! include 'COMMON.IOUNITS'
1782 ! include 'COMMON.HEADER'
1783 ! include 'COMMON.SBRIDGE'
1784 ! include 'COMMON.TEMPFAC'
1785 character(len=50) :: tytul
1786 character(len=1) :: chainid(10)=(/'A','B','C','D','E','F',&
1788 integer :: ica(nres)
1789 real(kind=8) :: etot,rmsd
1790 integer :: iatom,ichain,ires,i,j,iti,mnum,mnum1,boxxxshift(3)
1791 real(kind=8) :: boxxxx(3)
1792 write (ipdb,'(3a,1pe15.5,a,0pf7.2)') 'REMARK ',tytul(:20),&
1793 ' ENERGY ',etot,' RMS ',rmsd
1797 boxxxshift(1)=int(c(1,nnt)/boxxsize)
1798 ! if (c(1,nnt).lt.0) boxxxshift(1)=boxxxshift(1)-1
1799 boxxxshift(2)=int(c(2,nnt)/boxzsize)
1800 ! if (c(2,nnt).lt.0) boxxxshift(2)=boxxxshift(2)-1
1801 boxxxshift(3)=int(c(3,nnt)/boxzsize)
1802 ! if (c(3,nnt).lt.0) boxxxshift(3)=boxxxshift(3)-1
1811 if (iti.eq.ntyp1_molec(mnum)) then
1814 write (ipdb,'(a)') 'TER'
1821 if ((c(j,i).gt.0).and.(c(j,nnt).lt.0)) then
1822 c(j,i)=dmod(c(j,i),boxxxx(j))+(boxxxshift(j)-1)*boxxxx(j)
1823 else if ((c(j,i).lt.0).and.(c(j,nnt).gt.0)) then
1824 c(j,i)=dmod(c(j,i),boxxxx(j))+(boxxxshift(j)+1)*boxxxx(j)
1826 c(j,i)=dmod(c(j,i),boxxxx(j))+(boxxxshift(j))*boxxxx(j)
1830 write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),&
1831 ires,(c(j,i),j=1,3),1.0d0,tempfac(1,i)
1832 if ((iti.ne.10).and.(mnum.ne.5)) then
1834 write (ipdb,20) iatom,restyp(iti,mnum),chainid(ichain),&
1835 ires,(c(j,nres+i),j=1,3),1.0d0,tempfac(2,i)
1839 write (ipdb,'(a)') 'TER'
1843 if ((itype(i,mnum).eq.ntyp1).or.(mnum.eq.5)) cycle
1844 if (itype(i,mnum).eq.10 .and. itype(i+1,mnum1).ne.ntyp1_molec(mnum1)) then
1845 write (ipdb,30) ica(i),ica(i+1)
1846 else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum1).ne.ntyp1_molec(mnum1)) then
1847 write (ipdb,30) ica(i),ica(i+1),ica(i)+1
1848 else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) then
1849 write (ipdb,30) ica(i),ica(i)+1
1852 if ((itype(nct,molnum(nct)).ne.10).and.(molnum(nct).ne.5)) then
1853 write (ipdb,30) ica(nct),ica(nct)+1
1856 write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
1858 write (ipdb,'(a6)') 'ENDMDL'
1859 10 FORMAT ('ATOM',I7,' CA ',A3,1X,A1,I4,4X,3F8.3,2f6.2)
1860 20 FORMAT ('ATOM',I7,' CB ',A3,1X,A1,I4,4X,3F8.3,2f6.2)
1861 30 FORMAT ('CONECT',8I5)
1863 end subroutine pdboutC
1864 !-----------------------------------------------------------------------------
1865 subroutine cartout(igr,i,etot,free,rmsd,plik)
1866 ! implicit real*8 (a-h,o-z)
1867 ! include 'DIMENSIONS'
1868 ! include 'sizesclu.dat'
1869 ! include 'COMMON.IOUNITS'
1870 ! include 'COMMON.CHAIN'
1871 ! include 'COMMON.VAR'
1872 ! include 'COMMON.LOCAL'
1873 ! include 'COMMON.INTERACT'
1874 ! include 'COMMON.NAMES'
1875 ! include 'COMMON.GEO'
1876 ! include 'COMMON.CLUSTER'
1877 integer :: igr,i,j,k
1878 real(kind=8) :: etot,free,rmsd
1879 character(len=80) :: plik
1880 open (igeom,file=plik,position='append')
1881 write (igeom,'(2e15.5,f10.5,$)') etot,free,rmsd
1882 write (igeom,'(i4,$)') &
1883 nss_all(i),(ihpb_all(j,i),jhpb_all(j,i),j=1,nss_all(i))
1884 write (igeom,'(i10)') iscore(i)
1885 write (igeom,'(8f10.5)') &
1886 ((allcart(k,j,i),k=1,3),j=1,nres),&
1887 ((allcart(k,j+nres,i),k=1,3),j=nnt,nct)
1889 end subroutine cartout
1890 !------------------------------------------------------------------------------
1891 ! subroutine alloc_clust_arrays(n_conf)
1896 ! allocate(diss(maxdist)) !(maxdist)
1897 !el allocate(energy(0:maxconf),totfree(0:maxconf)) !(0:maxconf)
1898 ! allocatable :: enetb !(max_ene,maxstr_proc)
1899 !el allocate(entfac(maxconf)) !(maxconf)
1900 ! allocatable :: totfree_gr !(maxgr)
1901 !el allocate(rcutoff(max_cut+1)) !(max_cut+1)
1903 ! allocatable :: licz,iass !(maxgr)
1904 ! allocatable :: nconf !(maxgr,maxingr)
1905 ! allocatable :: iass_tot !(maxgr,max_cut)
1906 ! allocatable :: list_conf !(maxconf)
1908 !el allocatable :: allcart !(3,maxres2,maxstr_proc)
1909 !el allocate(rmstb(maxconf)) !(maxconf)
1910 !el allocate(mult(nres)) !(maxres)
1911 !el allocatable :: nss_all !(maxstr_proc)
1912 !el allocatable :: ihpb_all,jhpb_all !(maxss,maxstr_proc)
1913 ! allocate(icc(n_conf),iscore(n_conf)) !(maxconf)
1916 ! allocatable :: tempfac !(2,maxres)
1919 !el allocate(beta_h(maxT)) !(maxT)
1921 ! end subroutine alloc_clust_arrays
1922 !-----------------------------------------------------------------------------
1923 !-----------------------------------------------------------------------------