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)
569 if (from_cart .and. .not. from_bx .and. .not. from_cx) then
571 read (intin,*,end=13,err=11) energy(icon),totfree(icon),&
573 nss_all(icon),(ihpb_all(ii,icon),jhpb_all(i,icon),&
574 i=1,nss_all(icon)),iscore(icon)
576 read (intin,*,end=13,err=11) energy(icon),rmstb(icon),&
577 nss_all(icon),(ihpb_all(ii,icon),jhpb_all(i,icon),&
578 i=1,nss_all(icon)),iscore(icon)
580 read (intin,'(8f10.5)',end=13,err=10) &
581 ((allcart(j,i,icon),j=1,3),i=1,nres),&
582 ((allcart(j,i+nres,icon),j=1,3),i=nnt,nct)
583 print *,icon,energy(icon),nss_all(icon),rmstb(icon)
585 read(intin,'(a80)',end=13,err=12) lineh
586 read(lineh(:5),*,err=8) ic
588 read(lineh(6:),*,err=8) energy(icon)
590 read(lineh(6:),*,err=8) energy(icon)
594 print *,'error, assuming e=1d10',lineh
598 !old read(lineh(18:),*,end=13,err=11) nss_all(icon)
599 ii = index(lineh(15:)," ")+15
600 read(lineh(ii:),*,end=13,err=11) nss_all(icon)
601 IF (NSS_all(icon).LT.9) THEN
602 read (lineh(20:),*,end=102) &
603 (IHPB_all(I,icon),JHPB_all(I,icon),I=1,NSS_all(icon)),&
606 read (lineh(20:),*,end=102) &
607 (IHPB_all(I,icon),JHPB_all(I,icon),I=1,8)
608 read (intin,*) (IHPB_all(I,icon),JHPB_all(I,icon),&
609 I=9,NSS_all(icon)),iscore(icon)
614 PRINT *,'IC:',IC,' ENERGY:',ENERGY(ICON)
615 call read_angles(intin,*13)
617 phiall(i,icon)=phi(i)
618 thetall(i,icon)=theta(i)
619 alphall(i,icon)=alph(i)
620 omall(i,icon)=omeg(i)
626 ! CALCULATE DISTANCES
628 10 print *,'something wrong with angles'
630 11 print *,'something wrong with NSS',nss
632 12 print *,'something wrong with header'
639 open (icbase,file=bprotfiles,status="unknown",&
640 form="unformatted",access="direct",recl=lenrec)
641 ! Read conformations from binary DA files (one per batch) and write them to
642 ! a binary DA scratchfile.
646 write (liczba,'(bz,i3.3)') me
647 IF (ME.EQ.MASTER) THEN
648 ! Only the master reads the database; it'll send it to the other procs
656 open (intin,file=intinname,status="old",form="unformatted",&
657 access="direct",recl=lenrec_in)
659 else if (from_cx) then
660 #if (defined(AIX) && !defined(JUBL))
661 call xdrfopen_(ixdrf,intinname, "r", iret)
663 call xdrfopen(ixdrf,intinname, "r", iret)
666 write (iout,*) "xdrfopen: iret",iret
668 write (iout,*) "Error: coordinate file ",&
669 intinname(:ilen(intinname))," does not exist."
672 call MPI_ABORT(MPI_COMM_WORLD,IERROR,ERRCODE)
677 write (iout,*) "Error: coordinate format not specified"
680 call MPI_ABORT(MPI_COMM_WORLD,IERROR,ERRCODE)
688 write (iout,*) "Opening file ",intinname(:ilen(intinname))
689 write (iout,*) "lenrec",lenrec_in
693 ! write (iout,*) "maxconf",maxconf
697 !el if (i.gt.maxconf) then
698 !el write (iout,*) "Error: too many conformations ",&
699 !el "(",maxconf,") maximum."
701 !el call MPI_Abort(MPI_COMM_WORLD,errcode,ierror)
705 ! write (iout,*) "i",i
708 read(intin,err=101,end=101) &
709 ((csingle(l,k),l=1,3),k=1,nres),&
710 ((csingle(l,k+nres),l=1,3),k=nnt,nct),&
711 nss,(ihpb(k),jhpb(k),k=1,nss),&
713 entfac(jj+1),rmstb(jj+1),iscor
720 #if (defined(AIX) && !defined(JUBL))
721 call xdrf3dfcoord_(ixdrf, csingle, itmp, prec, iret)
722 if (iret.eq.0) goto 101
723 call xdrfint_(ixdrf, nss, iret)
724 if (iret.eq.0) goto 101
726 call xdrfint_(ixdrf, ihpb(j), iret)
727 if (iret.eq.0) goto 101
728 call xdrfint_(ixdrf, jhpb(j), iret)
729 if (iret.eq.0) goto 101
731 call xdrffloat_(ixdrf,reini,iret)
732 if (iret.eq.0) goto 101
733 call xdrffloat_(ixdrf,refree,iret)
734 if (iret.eq.0) goto 101
735 call xdrffloat_(ixdrf,rmsdev,iret)
736 if (iret.eq.0) goto 101
737 call xdrfint_(ixdrf,iscor,iret)
738 if (iret.eq.0) goto 101
740 ! write (iout,*) "calling xdrf3dfcoord"
741 call xdrf3dfcoord(ixdrf, csingle, itmp, prec, iret)
742 ! write (iout,*) "iret",iret
744 if (iret.eq.0) goto 101
745 call xdrfint(ixdrf, nss, iret)
746 write (iout,*) "iret",iret
747 ! write (iout,*) "nss",nss,i,"TUTU"
749 if (iret.eq.0) goto 101
751 call xdrfint(ixdrf, ihpb(k), iret)
752 if (iret.eq.0) goto 101
753 call xdrfint(ixdrf, jhpb(k), iret)
754 if (iret.eq.0) goto 101
755 ! write(iout,*), ihpb(k),jhpb(k),"TUTU"
757 call xdrffloat(ixdrf,reini,iret)
758 if (iret.eq.0) goto 101
759 ! write(iout,*) "TUTU", reini
760 call xdrffloat(ixdrf,refree,iret)
761 ! write(iout,*) "TUTU", refree
762 if (iret.eq.0) goto 101
763 call xdrffloat(ixdrf,rmsdev,iret)
764 if (iret.eq.0) goto 101
765 call xdrfint(ixdrf,iscor,iret)
766 if (iret.eq.0) goto 101
778 c(l,nres+k)=csingle(l,nres+k-nnt+1)
783 write (iout,'(5hREAD ,i5,3f15.4,i10)') &
784 jj+1,energy(jj+1),entfac(jj+1),&
786 write (iout,*) "Conformation",jjj+1,jj+1
787 write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
788 write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
791 call add_new_cconf(jjj,jj,jj_old,icount,Next)
794 write (iout,*) i-1," conformations read from DA file ",&
795 intinname(:ilen(intinname))
796 write (iout,*) jj," conformations read so far"
800 #if (defined(AIX) && !defined(JUBL))
801 call xdrfclose_(ixdrf, iret)
803 call xdrfclose(ixdrf, iret)
808 write (iout,*) "jj_old",jj_old," jj",jj
810 call write_and_send_cconf(icount,jj_old,jj,Next)
811 call MPI_Send(0,1,MPI_INTEGER,Next,570,&
812 MPI_COMM_WORLD,IERROR)
815 call write_and_send_cconf(icount,jj_old,jj,Next)
817 t_acq = tcpu() - t_acq
819 write (iout,*) "Processor",me,&
820 " time for conformation read/send",t_acq
822 ! A worker gets the confs from the master and sends them to its neighbor
824 call receive_and_pass_cconf(icount,jj_old,jj,&
826 t_acq = tcpu() - t_acq
833 write(iout,*)"A total of",ncon," conformations read."
835 allocate(enetb(1:max_ene,ncon)) !(max_ene,maxstr_proc)
837 ! Check if everyone has the same number of conformations
838 call MPI_Allgather(ncon,1,MPI_INTEGER,&
839 ntot_all(0),1,MPI_INTEGER,MPI_Comm_World,IERROR)
843 if (ncon.ne.ntot_all(i)) then
844 write (iout,*) "Number of conformations at processor",i,&
845 " differs from that at processor",me,&
853 write (iout,*) "Number of conformations read by processors"
856 write (iout,'(8i10)') i,ntot_all(i)
858 write (iout,*) "Calculation terminated."
864 1111 write(iout,*) "Error opening coordinate file ",&
865 intinname(:ilen(intinname))
868 end subroutine read_coords
869 !------------------------------------------------------------------------------
870 subroutine add_new_cconf(jjj,jj,jj_old,icount,Next)
872 use geometry_data, only: vbld,rad2deg,theta,phi,alph,omeg,deg2rad
873 use energy_data, only: itel,itype,dsc,max_ene,molnum
874 use control_data, only: symetr
875 use geometry, only: int_from_cart1
877 ! include "DIMENSIONS"
878 ! include "sizesclu.dat"
879 ! include "COMMON.CLUSTER"
880 ! include "COMMON.CONTROL"
881 ! include "COMMON.CHAIN"
882 ! include "COMMON.INTERACT"
883 ! include "COMMON.LOCAL"
884 ! include "COMMON.IOUNITS"
885 ! include "COMMON.NAMES"
886 ! include "COMMON.VAR"
887 ! include "COMMON.SBRIDGE"
888 ! include "COMMON.GEO"
889 integer :: i,j,jj,jjj,jj_old,icount,k,kk,l,ii,ib,&
890 nn,nn1,inan,Next,itj,chalen,mnum
891 real(kind=8) :: etot,energia(0:max_ene)
893 chalen=int((nct-nnt+2)/symetr)
894 call int_from_cart1(.false.)
897 write (iout,*) "Check atom",j,itype(j,mnum),vbld(j)
899 if (itype(j,mnum).eq.ntyp1_molec(mnum)) cycle
900 if (itype(j-1,molnum(j-1)).eq.ntyp1_molec(molnum(j-1))) cycle
902 if ((vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0).and.(mnum.ne.5)) then
904 if (itel(j).ne.0 .and. itel(j-1).ne.0) then
905 write (iout,*) "Conformation",jjj,jj+1,itype(j-1,molnum(j-1)),itype(j,mnum)
906 write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),itel(j),&
908 write (iout,*) "The Cartesian geometry is:"
909 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
910 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
911 write (iout,*) "The internal geometry is:"
912 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
913 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
914 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
915 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
916 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
917 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
919 "This conformation WILL NOT be added to the database."
929 if (itype(j,1).ne.10 .and. (itype(j,mnum).ne.ntyp1_molec(mnum)) &
930 .and. (vbld(nres+j)-dsc(iabs(itj))) &
932 write (iout,*) "Conformation",jjj,jj+1
933 write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j)
934 write (iout,*) "The Cartesian geometry is:"
935 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
936 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
937 write (iout,*) "The internal geometry is:"
938 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
939 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
940 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
941 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
942 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
943 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
945 "This conformation WILL NOT be added to the database."
950 if (theta(j).le.0.0d0) then
952 "Zero theta angle(s) in conformation",jjj,jj+1
953 write (iout,*) "The Cartesian geometry is:"
954 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
955 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
956 write (iout,*) "The internal geometry is:"
957 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
958 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
959 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
960 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
961 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
962 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
964 "This conformation WILL NOT be added to the database."
967 if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad
971 write (iout,*) "Conformation",jjj,jj
972 write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
973 write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
974 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
975 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
976 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
977 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
978 write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct)
979 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
980 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
981 write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
982 write (iout,'(e15.5,16i5)') entfac(icount+1)
983 ! & iscore(icount+1,0)
986 call store_cconf_from_file(jj,icount)
987 if (icount.eq.maxstr_proc) then
989 write (iout,* ) "jj_old",jj_old," jj",jj
991 call write_and_send_cconf(icount,jj_old,jj,Next)
996 end subroutine add_new_cconf
997 !------------------------------------------------------------------------------
998 subroutine store_cconf_from_file(jj,icount)
1000 use energy_data, only: ihpb,jhpb
1002 ! include "DIMENSIONS"
1003 ! include "sizesclu.dat"
1004 ! include "COMMON.CLUSTER"
1005 ! include "COMMON.CHAIN"
1006 ! include "COMMON.SBRIDGE"
1007 ! include "COMMON.INTERACT"
1008 ! include "COMMON.IOUNITS"
1009 ! include "COMMON.VAR"
1010 integer :: i,j,jj,icount
1011 ! Store the conformation that has been read in
1014 allcart(j,i,icount)=c(j,i)
1019 ihpb_all(i,icount)=ihpb(i)
1020 jhpb_all(i,icount)=jhpb(i)
1023 end subroutine store_cconf_from_file
1024 !------------------------------------------------------------------------------
1025 subroutine write_and_send_cconf(icount,jj_old,jj,Next)
1028 ! include "DIMENSIONS"
1029 ! include "sizesclu.dat"
1034 ! include "COMMON.MPI"
1036 ! include "COMMON.CHAIN"
1037 ! include "COMMON.SBRIDGE"
1038 ! include "COMMON.INTERACT"
1039 ! include "COMMON.IOUNITS"
1040 ! include "COMMON.CLUSTER"
1041 ! include "COMMON.VAR"
1042 integer :: icount,jj_old,jj,Next
1043 ! Write the structures to a scratch file
1045 ! Master sends the portion of conformations that have been read in to the neighbor
1047 write (iout,*) "Processor",me," entered WRITE_AND_SEND_CONF"
1050 call MPI_Send(icount,1,MPI_INTEGER,Next,570,MPI_COMM_WORLD,IERROR)
1051 call MPI_Send(nss_all(1),icount,MPI_INTEGER,&
1052 Next,571,MPI_COMM_WORLD,IERROR)
1053 call MPI_Send(ihpb_all(1,1),icount,MPI_INTEGER,&
1054 Next,572,MPI_COMM_WORLD,IERROR)
1055 call MPI_Send(jhpb_all(1,1),icount,MPI_INTEGER,&
1056 Next,573,MPI_COMM_WORLD,IERROR)
1057 call MPI_Send(rmstb(jj_old),icount,MPI_DOUBLE_PRECISION,&
1058 Next,577,MPI_COMM_WORLD,IERROR)
1059 call MPI_Send(entfac(jj_old),icount,MPI_DOUBLE_PRECISION,&
1060 Next,579,MPI_COMM_WORLD,IERROR)
1061 call MPI_Send(allcart(1,1,1),3*icount*2*nres,&
1062 MPI_REAL,Next,580,MPI_COMM_WORLD,IERROR)
1064 call dawrite_ccoords(jj_old,jj,icbase)
1066 end subroutine write_and_send_cconf
1067 !------------------------------------------------------------------------------
1069 subroutine receive_and_pass_cconf(icount,jj_old,jj,Previous,Next)
1073 ! include "DIMENSIONS"
1074 ! include "sizesclu.dat"
1076 integer :: IERROR !,STATUS(MPI_STATUS_SIZE)
1077 ! include "COMMON.MPI"
1078 ! include "COMMON.CHAIN"
1079 ! include "COMMON.SBRIDGE"
1080 ! include "COMMON.INTERACT"
1081 ! include "COMMON.IOUNITS"
1082 ! include "COMMON.VAR"
1083 ! include "COMMON.GEO"
1084 ! include "COMMON.CLUSTER"
1085 integer :: i,j,k,icount,jj_old,jj,Previous,Next
1088 write (iout,*) "Processor",me," entered RECEIVE_AND_PASS_CONF"
1091 do while (icount.gt.0)
1092 call MPI_Recv(icount,1,MPI_INTEGER,Previous,570,MPI_COMM_WORLD,&
1094 call MPI_Send(icount,1,MPI_INTEGER,Next,570,MPI_COMM_WORLD,&
1097 write (iout,*) "Processor",me," icount",icount
1099 if (icount.eq.0) return
1100 call MPI_Recv(nss_all(1),icount,MPI_INTEGER,&
1101 Previous,571,MPI_COMM_WORLD,STATUS,IERROR)
1102 call MPI_Send(nss_all(1),icount,MPI_INTEGER,&
1103 Next,571,MPI_COMM_WORLD,IERROR)
1104 call MPI_Recv(ihpb_all(1,1),icount,MPI_INTEGER,&
1105 Previous,572,MPI_COMM_WORLD,STATUS,IERROR)
1106 call MPI_Send(ihpb_all(1,1),icount,MPI_INTEGER,&
1107 Next,572,MPI_COMM_WORLD,IERROR)
1108 call MPI_Recv(jhpb_all(1,1),icount,MPI_INTEGER,&
1109 Previous,573,MPI_COMM_WORLD,STATUS,IERROR)
1110 call MPI_Send(jhpb_all(1,1),icount,MPI_INTEGER,&
1111 Next,573,MPI_COMM_WORLD,IERROR)
1112 call MPI_Recv(rmstb(jj_old),icount,MPI_DOUBLE_PRECISION,&
1113 Previous,577,MPI_COMM_WORLD,STATUS,IERROR)
1114 call MPI_Send(rmstb(jj_old),icount,MPI_DOUBLE_PRECISION,&
1115 Next,577,MPI_COMM_WORLD,IERROR)
1116 call MPI_Recv(entfac(jj_old),icount,MPI_DOUBLE_PRECISION,&
1117 Previous,579,MPI_COMM_WORLD,STATUS,IERROR)
1118 call MPI_Send(entfac(jj_old),icount,MPI_DOUBLE_PRECISION,&
1119 Next,579,MPI_COMM_WORLD,IERROR)
1120 call MPI_Recv(allcart(1,1,1),3*icount*2*nres,&
1121 MPI_REAL,Previous,580,MPI_COMM_WORLD,STATUS,IERROR)
1122 call MPI_Send(allcart(1,1,1),3*icount*2*nres,&
1123 MPI_REAL,Next,580,MPI_COMM_WORLD,IERROR)
1125 call dawrite_ccoords(jj_old,jj,icbase)
1128 write (iout,*) "Processor",me," received",icount," conformations"
1130 write (iout,'(8f10.4)') (allcart(l,k,i),l=1,3,k=1,nres)
1131 write (iout,'(8f10.4)')((allcart(l,k,i+nres),l=1,3,k=nnt,nct)
1132 write (iout,'(e15.5,16i5)') entfac(i)
1137 end subroutine receive_and_pass_cconf
1139 !------------------------------------------------------------------------------
1140 subroutine daread_ccoords(istart_conf,iend_conf)
1142 use energy_data, only: dyn_ss
1145 ! include "DIMENSIONS"
1146 ! include "sizesclu.dat"
1150 ! include "COMMON.MPI"
1152 ! include "COMMON.CHAIN"
1153 ! include "COMMON.CLUSTER"
1154 ! include "COMMON.IOUNITS"
1155 ! include "COMMON.INTERACT"
1156 ! include "COMMON.VAR"
1157 ! include "COMMON.SBRIDGE"
1158 ! include "COMMON.GEO"
1159 integer :: istart_conf,iend_conf
1160 integer :: i,j,ij,ii,iii
1162 character(len=16) :: form,acc
1163 character(len=32) :: nam
1165 ! Read conformations off a DA scratchfile.
1168 write (iout,*) "DAREAD_COORDS"
1169 write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
1170 inquire(unit=icbase,name=nam,recl=len,form=form,access=acc)
1171 write (iout,*) "len=",len," form=",form," acc=",acc
1172 write (iout,*) "nam=",nam
1175 do ii=istart_conf,iend_conf
1176 ij = ii - istart_conf + 1
1179 write (iout,*) "Reading binary file, record",iii," ii",ii
1183 read(icbase,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres), &
1184 ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
1185 entfac(ii),rmstb(ii)
1187 read(icbase,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres),&
1188 ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
1189 nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss),&
1190 entfac(ii),rmstb(ii)
1194 write (iout,*) ii,iii,ij,entfac(ii)
1195 write (iout,'(8f10.5)') ((allcart(j,i,ij),j=1,3),i=1,nres)
1196 write (iout,'(8f10.4)') ((allcart(j,i,ij),j=1,3),&
1197 i=nnt+nres,nct+nres)
1198 write (iout,'(2e15.5)') entfac(ij)
1199 write (iout,'(16i5)') nss_all(ij),(ihpb_all(i,ij),&
1200 jhpb_all(i,ij),i=1,nss)
1205 end subroutine daread_ccoords
1206 !------------------------------------------------------------------------------
1207 subroutine dawrite_ccoords(istart_conf,iend_conf,unit_out)
1210 ! include "DIMENSIONS"
1211 ! include "sizesclu.dat"
1212 use energy_data, only: dyn_ss
1217 ! include "COMMON.MPI"
1219 ! include "COMMON.CHAIN"
1220 ! include "COMMON.INTERACT"
1221 ! include "COMMON.IOUNITS"
1222 ! include "COMMON.VAR"
1223 ! include "COMMON.SBRIDGE"
1224 ! include "COMMON.GEO"
1225 ! include "COMMON.CLUSTER"
1226 integer :: istart_conf,iend_conf
1227 integer :: i,j,ii,ij,iii,unit_out
1229 character(len=16) :: form,acc
1230 character(len=32) :: nam
1232 ! Write conformations to a DA scratchfile.
1235 write (iout,*) "DAWRITE_COORDS"
1236 write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
1237 write (iout,*) "lenrec",lenrec
1238 inquire(unit=unit_out,name=nam,recl=len,form=form,access=acc)
1239 write (iout,*) "len=",len," form=",form," acc=",acc
1240 write (iout,*) "nam=",nam
1243 do ii=istart_conf,iend_conf
1245 ij = ii - istart_conf + 1
1247 write (iout,*) "Writing binary file, record",iii," ii",ii
1251 write(unit_out,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres),&
1252 ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
1253 entfac(ii),rmstb(ii)
1255 write(unit_out,rec=iii) ((allcart(j,i,ij),j=1,3),i=1,nres), &
1256 ((allcart(j,i,ij),j=1,3),i=nnt+nres,nct+nres),&
1257 nss_all(ij),(ihpb_all(i,ij),jhpb_all(i,ij),i=1,nss_all(ij)),&
1258 entfac(ii),rmstb(ii)
1262 write (iout,'(8f10.5)') ((allcart(j,i,ij),j=1,3),i=1,nres)
1263 write (iout,'(8f10.4)') ((allcart(j,i,ij),j=1,3),i=nnt+nres,&
1265 write (iout,'(2e15.5)') entfac(ij)
1266 write (iout,'(16i5)') nss_all(ij),(ihpb(i,ij),jhpb(i,ij),i=1,&
1272 end subroutine dawrite_ccoords
1273 !-----------------------------------------------------------------------------
1275 !-----------------------------------------------------------------------------
1276 subroutine read_control
1278 ! Read molecular data
1280 use energy_data, only: rescale_mode,distchainmax,ipot !,temp0
1281 use control_data, only: titel,outpdb,outmol2,refstr,pdbref,&
1282 iscode,symetr,punch_dist,print_dist,nstart,nend,&
1283 caonly,iopt,efree,lprint_cart,lprint_int,rlamb_ele,&
1284 r_cut_ele,nclust,tor_mode,scelemode,r_cut_mart,r_cut_ang,&
1286 use geometry_data, only:bordliptop,bordlipbot,&
1287 bufliptop,buflipbot,lipthick,lipbufthick
1289 ! include 'DIMENSIONS'
1290 ! include 'sizesclu.dat'
1291 ! include 'COMMON.IOUNITS'
1292 ! include 'COMMON.TIME1'
1293 ! include 'COMMON.SBRIDGE'
1294 ! include 'COMMON.CONTROL'
1295 ! include 'COMMON.CLUSTER'
1296 ! include 'COMMON.CHAIN'
1297 ! include 'COMMON.HEADER'
1298 ! include 'COMMON.FFIELD'
1299 ! include 'COMMON.FREE'
1300 ! include 'COMMON.INTERACT'
1301 character(len=320) :: controlcard !,ucase
1303 ! include 'COMMON.INFO'
1307 read (INP,'(a80)') titel
1308 call card_concat(controlcard,.true.)
1310 call readi(controlcard,'NRES',nres_molec(1),0)
1311 call readi(controlcard,'NRES_NUCL',nres_molec(2),0)
1312 call readi(controlcard,'NRES_CAT',nres_molec(5),0)
1315 nres=nres_molec(i)+nres
1317 print *,"TU",nres_molec(:)
1319 ! call alloc_clust_arrays
1320 allocate(rcutoff(max_cut+1)) !(max_cut+1)
1321 allocate(beta_h(maxT)) !(maxT)
1322 allocate(mult(nres)) !(maxres)
1325 call readi(controlcard,'RESCALE',rescale_mode,2)
1326 call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
1327 write (iout,*) "DISTCHAINMAX",distchainmax
1328 call readi(controlcard,'PDBOUT',outpdb,0)
1329 call readi(controlcard,'MOL2OUT',outmol2,0)
1330 refstr=(index(controlcard,'REFSTR').gt.0)
1331 write (iout,*) "REFSTR",refstr
1332 pdbref=(index(controlcard,'PDBREF').gt.0)
1333 iscode=index(controlcard,'ONE_LETTER')
1334 tree=(index(controlcard,'MAKE_TREE').gt.0)
1335 min_var=(index(controlcard,'MINVAR').gt.0)
1336 plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
1337 punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
1338 call readi(controlcard,'NCUT',ncut,1)
1339 call readi(controlcard,'SYM',symetr,1)
1340 write (iout,*) 'sym', symetr
1341 call readi(controlcard,'NSTART',nstart,0)
1342 call reada(controlcard,'BOXX',boxxsize,100.0d0)
1343 call reada(controlcard,'BOXY',boxysize,100.0d0)
1344 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
1345 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
1346 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
1347 if (lipthick.gt.0.0d0) then
1348 bordliptop=(boxzsize+lipthick)/2.0
1349 bordlipbot=bordliptop-lipthick
1350 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
1351 write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
1352 buflipbot=bordlipbot+lipbufthick
1353 bufliptop=bordliptop-lipbufthick
1354 if ((lipbufthick*2.0d0).gt.lipthick) &
1355 write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
1356 endif !lipthick.gt.0
1357 write(iout,*) "bordliptop=",bordliptop
1358 write(iout,*) "bordlipbot=",bordlipbot
1359 write(iout,*) "bufliptop=",bufliptop
1360 write(iout,*) "buflipbot=",buflipbot
1362 call readi(controlcard,'NCLUST',nclust,5)
1363 ! ions=index(controlcard,"IONS").gt.0
1364 call readi(controlcard,"SCELEMODE",scelemode,0)
1365 print *,"SCELE",scelemode
1366 call readi(controlcard,'TORMODE',tor_mode,0)
1367 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1368 write(iout,*) "torsional and valence angle mode",tor_mode
1369 call reada(controlcard,"R_CUT_ELE",r_cut_ele,25.0d0)
1370 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
1371 write(iout,*) "R_CUT_ELE=",r_cut_ele
1372 call reada(controlcard,"R_CUT_MART",r_cut_mart,15.0d0)
1373 call reada(controlcard,"LAMBDA_MART",rlamb_mart,0.3d0)
1374 call reada(controlcard,"R_CUT_ANG",r_cut_ang,4.2d0)
1375 call readi(controlcard,'NEND',nend,0)
1376 call reada(controlcard,'ECUT',ecut,10.0d0)
1377 call reada(controlcard,'PROB',prob_limit,0.99d0)
1378 write (iout,*) "Probability limit",prob_limit
1379 lgrp=(index(controlcard,'LGRP').gt.0)
1380 caonly=(index(controlcard,'CA_ONLY').gt.0)
1381 print_dist=(index(controlcard,'PRINT_DIST').gt.0)
1382 call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
1383 call readi(controlcard,'IOPT',iopt,2)
1384 lside = index(controlcard,"SIDE").gt.0
1385 efree = index(controlcard,"EFREE").gt.0
1386 call readi(controlcard,'NTEMP',nT,1)
1387 write (iout,*) "nT",nT
1388 !el call reada(controlcard,'TEMP0',temp0,300.0d0) !el
1389 call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
1390 write (iout,*) "nT",nT
1391 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
1393 beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
1395 write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
1396 lprint_cart=index(controlcard,"PRINT_CART") .gt.0
1397 lprint_int=index(controlcard,"PRINT_INT") .gt.0
1400 end subroutine read_control
1401 !-----------------------------------------------------------------------------
1404 ! Read molecular data.
1406 use geometry_data, only: nsup,cref,nres0,nstart_sup,nstart_seq,dc
1407 use energy_data!, only: wsc,wscp,welec,wbond,wstrain,wtor,wtor_d,&
1408 ! wang,wscloc,wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,&
1409 ! wturn3,wturn4,wturn6,wvdwpp,weights
1410 use control_data, only: titel,nstart,nend,pdbref,refstr,iscode,&
1412 use geometry, only: chainbuild,alloc_geo_arrays
1413 use energy, only: alloc_ener_arrays
1414 use io_base, only: rescode
1415 use control, only: setup_var,init_int_table
1416 use conform_compar, only: contact
1418 ! include 'DIMENSIONS'
1419 ! include 'COMMON.IOUNITS'
1420 ! include 'COMMON.GEO'
1421 ! include 'COMMON.VAR'
1422 ! include 'COMMON.INTERACT'
1423 ! include 'COMMON.LOCAL'
1424 ! include 'COMMON.NAMES'
1425 ! include 'COMMON.CHAIN'
1426 ! include 'COMMON.FFIELD'
1427 ! include 'COMMON.SBRIDGE'
1428 ! include 'COMMON.HEADER'
1429 ! include 'COMMON.CONTROL'
1430 ! include 'COMMON.CONTACTS'
1431 ! include 'COMMON.TIME1'
1433 ! include 'COMMON.INFO'
1435 character(len=4) :: sequence(nres) !(maxres)
1436 character(len=800) :: weightcard
1438 real(kind=8) :: x(6*nres) !(maxvar)
1439 real(kind=8) :: ssscale
1440 integer :: itype_pdb(nres) !(maxres)
1442 integer :: i,j,kkk,mnum
1446 !el allocate(weights(n_ene))
1447 allocate(weights(max_ene))
1448 call alloc_geo_arrays
1449 call alloc_ener_arrays
1450 !-----------------------------
1451 allocate(c(3,2*nres+2)) !(3,maxres2+2) maxres2=2*maxres
1452 allocate(dc(3,0:2*nres+2)) !(3,0:maxres2)
1453 allocate(itype(nres+2,5)) !(maxres)
1454 allocate(molnum(nres+1))
1455 allocate(itel(nres+2))
1469 !--------------------------
1470 ! Read weights of the subsequent energy terms.
1471 call card_concat(weightcard,.true.)
1472 call reada(weightcard,'WSC',wsc,1.0d0)
1473 call reada(weightcard,'WLONG',wsc,wsc)
1474 call reada(weightcard,'WSCP',wscp,1.0d0)
1475 call reada(weightcard,'WELEC',welec,1.0D0)
1476 call reada(weightcard,'WVDWPP',wvdwpp,welec)
1477 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
1478 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
1479 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
1480 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
1481 call reada(weightcard,'WTURN3',wturn3,1.0D0)
1482 call reada(weightcard,'WTURN4',wturn4,1.0D0)
1483 call reada(weightcard,'WTURN6',wturn6,1.0D0)
1484 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
1485 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
1486 call reada(weightcard,'WBOND',wbond,1.0D0)
1487 call reada(weightcard,'WTOR',wtor,1.0D0)
1488 call reada(weightcard,'WTORD',wtor_d,1.0D0)
1489 call reada(weightcard,'WANG',wang,1.0D0)
1490 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
1491 call reada(weightcard,'SCAL14',scal14,0.4D0)
1492 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
1493 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
1494 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
1495 call reada(weightcard,'WCATCAT',wcatcat,0.0d0)
1496 call reada(weightcard,'WCATPROT',wcatprot,0.0d0)
1497 call reada(weightcard,'TEMP0',temp0,300.0d0) !!! el
1498 if (index(weightcard,'SOFT').gt.0) ipot=6
1499 call reada(weightcard,'WVDWPP_NUCL',wvdwpp_nucl,0.0D0)
1500 call reada(weightcard,'WELPP',welpp,0.0d0)
1501 call reada(weightcard,'WVDWPSB',wvdwpsb,0.0d0)
1502 call reada(weightcard,'WELPSB',welpsb,0.0D0)
1503 call reada(weightcard,'WVDWSB',wvdwsb,0.0d0)
1504 call reada(weightcard,'WELSB',welsb,0.0D0)
1505 call reada(weightcard,'WBOND_NUCL',wbond_nucl,0.0D0)
1506 call reada(weightcard,'WANG_NUCL',wang_nucl,0.0D0)
1507 call reada(weightcard,'WSBLOC',wsbloc,0.0D0)
1508 call reada(weightcard,'WTOR_NUCL',wtor_nucl,0.0D0)
1509 ! print *,"KUR..",wtor_nucl
1510 call reada(weightcard,'WTORD_NUCL',wtor_d_nucl,0.0D0)
1511 call reada(weightcard,'WCORR_NUCL',wcorr_nucl,0.0D0)
1512 call reada(weightcard,'WCORR3_NUC',wcorr3_nucl,0.0D0)
1513 ! 12/1/95 Added weight for the multi-body term WCORR
1514 call reada(weightcard,'WCORRH',wcorr,1.0D0)
1515 call reada(weightcard,'WSCBASE',wscbase,0.0D0)
1516 call reada(weightcard,'WPEPBASE',wpepbase,1.0d0)
1517 call reada(weightcard,'WSCPHO',wscpho,0.0d0)
1518 call reada(weightcard,'WPEPPHO',wpeppho,0.0d0)
1519 call reada(weightcard,'WCATNUCL',wcatnucl,0.0d0)
1520 call reada(weightcard,'WCATTRAN',wcat_tran,0.0d0)
1521 call reada(weightcard,"D0CM",d0cm,3.78d0)
1522 call reada(weightcard,"AKCM",akcm,15.1d0)
1523 call reada(weightcard,"AKTH",akth,11.0d0)
1524 call reada(weightcard,"AKCT",akct,12.0d0)
1525 call reada(weightcard,"V1SS",v1ss,-1.08d0)
1526 call reada(weightcard,"V2SS",v2ss,7.61d0)
1527 call reada(weightcard,"V3SS",v3ss,13.7d0)
1528 call reada(weightcard,"EBR",ebr,-5.50D0)
1529 call reada(weightcard,"ATRISS",atriss,0.301D0)
1530 call reada(weightcard,"BTRISS",btriss,0.021D0)
1531 call reada(weightcard,"CTRISS",ctriss,1.001D0)
1532 call reada(weightcard,"DTRISS",dtriss,1.001D0)
1533 call reada(weightcard,"SSSCALE",ssscale,1.0D0)
1534 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
1535 call reada(weightcard,"LIPSCALE",lipscale,1.0D0)
1536 call reada(weightcard,"WTL",wliptran,1.0D0)
1538 call reada(weightcard,"HT",Ht,0.0D0)
1540 ss_depth=(ebr/wsc-0.25*eps(1,1))*ssscale
1541 Ht=(Ht/wsc-0.25*eps(1,1))*ssscale
1542 akcm=akcm/wsc*ssscale
1543 akth=akth/wsc*ssscale
1544 akct=akct/wsc*ssscale
1545 v1ss=v1ss/wsc*ssscale
1546 v2ss=v2ss/wsc*ssscale
1547 v3ss=v3ss/wsc*ssscale
1549 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
1552 if (wcorr4.gt.0.0d0) wcorr=wcorr4
1571 !el weights(19)=wsccor !!!!!!!!!!!!!!!!
1573 weights(26)=wvdwpp_nucl
1579 weights(32)=wbond_nucl
1580 weights(33)=wang_nucl
1582 weights(35)=wtor_nucl
1583 weights(36)=wtor_d_nucl
1584 weights(37)=wcorr_nucl
1585 weights(38)=wcorr3_nucl
1587 weights(42)=wcatprot
1592 weights(50)=wcatnucl
1593 weights(56)=wcat_tran
1594 weights(22)=wliptran
1596 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
1597 wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,&
1598 wturn4,wturn6,wsccor
1599 10 format (/'Energy-term weights (unscaled):'// &
1600 'WSCC= ',f10.6,' (SC-SC)'/ &
1601 'WSCP= ',f10.6,' (SC-p)'/ &
1602 'WELEC= ',f10.6,' (p-p electr)'/ &
1603 'WVDWPP= ',f10.6,' (p-p VDW)'/ &
1604 'WBOND= ',f10.6,' (stretching)'/ &
1605 'WANG= ',f10.6,' (bending)'/ &
1606 'WSCLOC= ',f10.6,' (SC local)'/ &
1607 'WTOR= ',f10.6,' (torsional)'/ &
1608 'WTORD= ',f10.6,' (double torsional)'/ &
1609 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ &
1610 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ &
1611 'WCORR4= ',f10.6,' (multi-body 4th order)'/ &
1612 'WCORR5= ',f10.6,' (multi-body 5th order)'/ &
1613 'WCORR6= ',f10.6,' (multi-body 6th order)'/ &
1614 'WTURN3= ',f10.6,' (turns, 3rd order)'/ &
1615 'WTURN4= ',f10.6,' (turns, 4th order)'/ &
1616 'WTURN6= ',f10.6,' (turns, 6th order)'/ &
1617 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
1619 if (wcorr4.gt.0.0d0) then
1620 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',&
1621 'between contact pairs of peptide groups'
1622 write (iout,'(2(a,f5.3/))') &
1623 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,&
1624 'Range of quenching the correlation terms:',2*delt_corr
1625 else if (wcorr.gt.0.0d0) then
1626 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',&
1627 'between contact pairs of peptide groups'
1629 write (iout,'(a,f8.3)') &
1630 'Scaling factor of 1,4 SC-p interactions:',scal14
1631 write (iout,'(a,f8.3)') &
1632 'General scaling factor of SC-p interactions:',scalscp
1633 r0_corr=cutoff_corr-delt_corr
1635 aad(i,1)=scalscp*aad(i,1)
1636 aad(i,2)=scalscp*aad(i,2)
1637 bad(i,1)=scalscp*bad(i,1)
1638 bad(i,2)=scalscp*bad(i,2)
1642 print *,'indpdb=',indpdb,' pdbref=',pdbref
1644 ! Read sequence if not taken from the pdb file.
1645 if (iscode.gt.0) then
1646 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
1648 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
1650 print *,nres_molec(:),nres
1651 ! Convert sequence to numeric code
1652 do i=1,nres_molec(1)
1655 itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
1657 do i=nres_molec(1)+1,nres_molec(1)+nres_molec(2)
1660 write (iout,*),i,sequence(i)
1661 itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
1664 do i=nres_molec(1)+nres_molec(2)+1,nres_molec(1)+nres_molec(2)+nres_molec(5)
1667 itype(i,mnum)=rescode(i,sequence(i),iscode,mnum)
1670 print '(20i4)',(itype(i,molnum(i)),i=1,nres)
1674 if (itype(i,molnum(i)).eq.ntyp1_molec(molnum(i))&
1675 .or. itype(i+1,molnum(i+1)).eq.ntyp1_molec(molnum(i+1))) then
1677 if (itype(i,molnum(i)).eq.ntyp1_molec(molnum(i))) then
1681 else if (iabs(itype(i+1,1)).ne.20) then
1683 else if (iabs(itype(i,1)).ne.20) then
1690 write (iout,*) "ITEL"
1692 write (iout,*) i,itype(i,molnum(i)),itel(i)
1695 print *,'Call Read_Bridge.'
1699 print *,'NNT=',NNT,' NCT=',NCT
1700 if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2
1701 if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1
1702 if (nstart.lt.nnt) nstart=nnt
1703 if (nend.gt.nct .or. nend.eq.0) nend=nct
1704 write (iout,*) "nstart",nstart," nend",nend
1707 ! read(inp,'(a)') pdbfile
1708 ! write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
1709 ! open(ipdbin,file=pdbfile,status='old',err=33)
1711 ! 33 write (iout,'(a)') 'Error opening PDB file.'
1714 ! print *,'Begin reading pdb data'
1716 ! print *,'Finished reading pdb data'
1717 ! write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
1719 ! itype_pdb(i)=itype(i)
1722 ! write (iout,'(a,i3)') 'nsup=',nsup
1724 ! if (nsup.le.(nct-nnt+1)) then
1725 ! do i=0,nct-nnt+1-nsup
1726 ! if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1731 ! write (iout,'(a)')
1732 ! & 'Error - sequences to be superposed do not match.'
1735 ! do i=0,nsup-(nct-nnt+1)
1736 ! if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1738 ! nstart_sup=nstart_sup+i
1743 ! write (iout,'(a)')
1744 ! & 'Error - sequences to be superposed do not match.'
1747 ! write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1748 ! & ' nstart_seq=',nstart_seq
1750 write(iout,*)"przed ini_int_tab"
1752 write(iout,*)"po ini_int_tab"
1753 write(iout,*)"przed setup var"
1755 write(iout,*)"po setup var"
1756 if (ns.gt.0.and.dyn_ss) then
1760 forcon(i-nss)=forcon(i)
1767 ! call hpb_partition
1769 dyn_ss_mask(iss(i))=.true.
1773 write (iout,*) "molread: REFSTR",refstr
1775 if (.not.pdbref) then
1776 call read_angles(inp,*38)
1778 38 write (iout,'(a)') 'Error reading reference structure.'
1780 call mp_stopall(Error_Msg)
1782 stop 'Error reading reference structure'
1791 cref(j,i,kkk)=c(j,i)
1795 call contact(.true.,ncont_ref,icont_ref)
1798 end subroutine molread
1799 !-----------------------------------------------------------------------------
1800 subroutine openunits
1802 ! include 'DIMENSIONS'
1803 use control_data, only: from_cx,from_bx,from_cart
1807 character(len=3) :: liczba
1808 ! include "COMMON.MPI"
1810 ! include 'COMMON.IOUNITS'
1811 ! include 'COMMON.CONTROL'
1812 integer :: lenpre,lenpot !,ilen
1814 character(len=16) :: cformat,cprint
1815 ! character(len=16) ucase
1816 integer :: lenint,lenout
1817 call getenv('INPUT',prefix)
1818 call getenv('OUTPUT',prefout)
1819 call getenv('INTIN',prefintin)
1820 call getenv('COORD',cformat)
1821 call getenv('PRINTCOOR',cprint)
1822 call getenv('SCRATCHDIR',scratchdir)
1825 if (index(ucase(cformat),'CX').gt.0) then
1831 lenout=ilen(prefout)
1832 lenint=ilen(prefintin)
1833 ! Get the names and open the input files
1834 open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
1836 write (liczba,'(bz,i3.3)') me
1837 outname=prefout(:lenout)//'_clust.out_'//liczba
1839 outname=prefout(:lenout)//'_clust.out'
1842 intinname=prefintin(:lenint)//'.bx'
1843 else if (from_cx) then
1844 intinname=prefintin(:lenint)//'.cx'
1846 intinname=prefintin(:lenint)//'.int'
1848 rmsname=prefintin(:lenint)//'.rms'
1849 open (jplot,file=prefout(:ilen(prefout))//'.tex',&
1851 open (jrms,file=rmsname,status='unknown')
1852 open(iout,file=outname,status='unknown')
1853 ! Get parameter filenames and open the parameter files.
1854 call getenv('BONDPAR',bondname)
1855 open (ibond,file=bondname,status='old')
1856 call getenv('THETPAR',thetname)
1857 open (ithep,file=thetname,status='old')
1858 call getenv('ROTPAR',rotname)
1859 open (irotam,file=rotname,status='old')
1860 call getenv('TORPAR',torname)
1861 open (itorp,file=torname,status='old')
1862 call getenv('TORDPAR',tordname)
1863 open (itordp,file=tordname,status='old')
1864 call getenv('FOURIER',fouriername)
1865 open (ifourier,file=fouriername,status='old')
1866 call getenv('ELEPAR',elename)
1867 open (ielep,file=elename,status='old')
1868 call getenv('SIDEPAR',sidename)
1869 open (isidep,file=sidename,status='old')
1870 call getenv('SIDEP',sidepname)
1871 open (isidep1,file=sidepname,status="old")
1872 call getenv('SCCORPAR',sccorname)
1873 open (isccor,file=sccorname,status="old")
1874 call getenv('BONDPAR_NUCL',bondname_nucl)
1875 open (ibond_nucl,file=bondname_nucl,status='old')
1876 call getenv('THETPAR_NUCL',thetname_nucl)
1877 open (ithep_nucl,file=thetname_nucl,status='old')
1878 call getenv('ROTPAR_NUCL',rotname_nucl)
1879 open (irotam_nucl,file=rotname_nucl,status='old')
1880 call getenv('TORPAR_NUCL',torname_nucl)
1881 open (itorp_nucl,file=torname_nucl,status='old')
1882 call getenv('TORDPAR_NUCL',tordname_nucl)
1883 open (itordp_nucl,file=tordname_nucl,status='old')
1884 call getenv('SIDEPAR_NUCL',sidename_nucl)
1885 open (isidep_nucl,file=sidename_nucl,status='old')
1886 call getenv('SIDEPAR_SCBASE',sidename_scbase)
1887 open (isidep_scbase,file=sidename_scbase,status='old')
1888 call getenv('PEPPAR_PEPBASE',pepname_pepbase)
1889 open (isidep_pepbase,file=pepname_pepbase,status='old')
1890 call getenv('SCPAR_PHOSPH',pepname_scpho)
1891 open (isidep_scpho,file=pepname_scpho,status='old')
1892 call getenv('PEPPAR_PHOSPH',pepname_peppho)
1893 open (isidep_peppho,file=pepname_peppho,status='old')
1896 call getenv('LIPTRANPAR',liptranname)
1897 open (iliptranpar,file=liptranname,status='old')
1898 call getenv('TUBEPAR',tubename)
1899 open (itube,file=tubename,status='old')
1900 call getenv('IONPAR',ionname)
1901 open (iion,file=ionname,status='old')
1902 call getenv('IONPAR_NUCL',ionnuclname)
1903 open (iionnucl,file=ionnuclname,status='old')
1904 call getenv('IONPAR_TRAN',iontranname)
1905 open (iiontran,file=iontranname,status='old')
1910 ! 8/9/01 In the newest version SCp interaction constants are read from a file
1911 ! Use -DOLDSCP to use hard-coded constants instead.
1913 call getenv('SCPPAR',scpname)
1914 open (iscpp,file=scpname,status='old')
1915 call getenv('SCPPAR_NUCL',scpname_nucl)
1916 open (iscpp_nucl,file=scpname_nucl,status='old')
1920 end subroutine openunits
1921 !-----------------------------------------------------------------------------
1923 !-----------------------------------------------------------------------------
1924 subroutine pdboutC(etot,rmsd,tytul)
1926 use energy_data, only: ihpb,jhpb,itype,molnum
1927 use energy, only:boxshift,to_box
1928 ! implicit real*8 (a-h,o-z)
1929 ! include 'DIMENSIONS'
1930 ! include 'COMMON.CONTROL'
1931 ! include 'COMMON.CHAIN'
1932 ! include 'COMMON.INTERACT'
1933 ! include 'COMMON.NAMES'
1934 ! include 'COMMON.IOUNITS'
1935 ! include 'COMMON.HEADER'
1936 ! include 'COMMON.SBRIDGE'
1937 ! include 'COMMON.TEMPFAC'
1938 character(len=50) :: tytul
1939 character(len=1) :: chainid(10)=(/'A','B','C','D','E','F',&
1941 integer :: ica(nres),k,iti1,ki
1942 real(kind=8) :: etot,rmsd
1943 integer :: iatom,ichain,ires,i,j,iti,mnum,mnum1,boxxxshift(3)
1944 real(kind=8) :: boxxxx(3),cbeg(3),Rdist(3)
1945 write (ipdb,'(3a,1pe15.5,a,0pf7.2)') 'REMARK ',tytul(:20),&
1946 ' ENERGY ',etot,' RMS ',rmsd
1950 ! boxxxshift(1)=int(c(1,nnt)/boxxsize)
1951 ! if (c(1,nnt).lt.0) boxxxshift(1)=boxxxshift(1)-1
1952 ! boxxxshift(2)=int(c(2,nnt)/boxzsize)
1953 ! if (c(2,nnt).lt.0) boxxxshift(2)=boxxxshift(2)-1
1954 ! boxxxshift(3)=int(c(3,nnt)/boxzsize)
1955 ! if (c(3,nnt).lt.0) boxxxshift(3)=boxxxshift(3)-1
1957 ! if (boxxxshift(i).gt.2) boxxxshift(i)=2
1958 ! if (boxxxshift(i).lt.-2) boxxxshift(i)=-2
1964 call to_box(cbeg(1),cbeg(2),cbeg(3))
1973 iti1=itype(i+1,mnum)
1975 if ((iti.eq.ntyp1_molec(mnum)).and.(iti1.eq.ntyp1_molec(mnum1))) cycle
1977 if (iti.eq.ntyp1_molec(mnum)) then
1980 write (ipdb,'(a)') 'TER'
1983 ! if (molnum(i).ge.1) then
1987 if (itype(i,mnum).ne.ntyp1_molec(mnum)) then
1993 if (itype(ki,mnum).eq.ntyp1_molec(mnum)) ki=ki+1
1994 if (itype(ki,mnum).eq.ntyp1_molec(mnum)) ki=ki+1
2000 call to_box(c(1,i),c(2,i),c(3,i))
2003 Rdist(k) = boxshift(c(k,i) - cbeg(k),boxxxx(k))
2004 c(k,i)=cbeg(k)+Rdist(k)
2006 if ((itype(i,molnum(i)).ne.10).and.(molnum(i).ne.5)) then
2007 call to_box(c(1,i+nres),c(2,i+nres),c(3,i+nres))
2009 Rdist(k) = boxshift(c(k,i+nres) - cbeg(k),boxxxx(k))
2010 c(k,i+nres)=cbeg(k)+Rdist(k)
2013 write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),&
2014 ires,(c(j,i),j=1,3),1.0d0,tempfac(1,i)
2015 if ((iti.ne.10).and.(mnum.ne.5)) then
2017 write (ipdb,20) iatom,restyp(iti,mnum),chainid(ichain),&
2018 ires,(c(j,nres+i),j=1,3),1.0d0,tempfac(2,i)
2022 write (ipdb,'(a)') 'TER'
2026 if ((itype(i,mnum).eq.ntyp1).or.(mnum.eq.5)) cycle
2027 if (itype(i,mnum).eq.10 .and. itype(i+1,mnum1).ne.ntyp1_molec(mnum1)) then
2028 write (ipdb,30) ica(i),ica(i+1)
2029 else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum1).ne.ntyp1_molec(mnum1)) then
2030 write (ipdb,30) ica(i),ica(i+1),ica(i)+1
2031 else if (itype(i,mnum).ne.10 .and. itype(i+1,mnum1).eq.ntyp1_molec(mnum1)) then
2032 write (ipdb,30) ica(i),ica(i)+1
2035 if ((itype(nct,molnum(nct)).ne.10).and.(molnum(nct).ne.5)) then
2036 write (ipdb,30) ica(nct),ica(nct)+1
2039 write (ipdb,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
2041 write (ipdb,'(a6)') 'ENDMDL'
2042 10 FORMAT ('ATOM',I7,' CA ',A3,1X,A1,I4,4X,3F8.3,2f6.2)
2043 20 FORMAT ('ATOM',I7,' CB ',A3,1X,A1,I4,4X,3F8.3,2f6.2)
2044 30 FORMAT ('CONECT',8I5)
2046 end subroutine pdboutC
2047 !-----------------------------------------------------------------------------
2048 subroutine cartout(igr,i,etot,free,rmsd,plik)
2049 ! implicit real*8 (a-h,o-z)
2050 ! include 'DIMENSIONS'
2051 ! include 'sizesclu.dat'
2052 ! include 'COMMON.IOUNITS'
2053 ! include 'COMMON.CHAIN'
2054 ! include 'COMMON.VAR'
2055 ! include 'COMMON.LOCAL'
2056 ! include 'COMMON.INTERACT'
2057 ! include 'COMMON.NAMES'
2058 ! include 'COMMON.GEO'
2059 ! include 'COMMON.CLUSTER'
2060 integer :: igr,i,j,k
2061 real(kind=8) :: etot,free,rmsd
2062 character(len=80) :: plik
2063 open (igeom,file=plik,position='append')
2064 write (igeom,'(2e15.5,f10.5,$)') etot,free,rmsd
2065 write (igeom,'(i4,$)') &
2066 nss_all(i),(ihpb_all(j,i),jhpb_all(j,i),j=1,nss_all(i))
2067 write (igeom,'(i10)') iscore(i)
2068 write (igeom,'(8f10.5)') &
2069 ((allcart(k,j,i),k=1,3),j=1,nres),&
2070 ((allcart(k,j+nres,i),k=1,3),j=nnt,nct)
2072 end subroutine cartout
2073 !------------------------------------------------------------------------------
2074 ! subroutine alloc_clust_arrays(n_conf)
2079 ! allocate(diss(maxdist)) !(maxdist)
2080 !el allocate(energy(0:maxconf),totfree(0:maxconf)) !(0:maxconf)
2081 ! allocatable :: enetb !(max_ene,maxstr_proc)
2082 !el allocate(entfac(maxconf)) !(maxconf)
2083 ! allocatable :: totfree_gr !(maxgr)
2084 !el allocate(rcutoff(max_cut+1)) !(max_cut+1)
2086 ! allocatable :: licz,iass !(maxgr)
2087 ! allocatable :: nconf !(maxgr,maxingr)
2088 ! allocatable :: iass_tot !(maxgr,max_cut)
2089 ! allocatable :: list_conf !(maxconf)
2091 !el allocatable :: allcart !(3,maxres2,maxstr_proc)
2092 !el allocate(rmstb(maxconf)) !(maxconf)
2093 !el allocate(mult(nres)) !(maxres)
2094 !el allocatable :: nss_all !(maxstr_proc)
2095 !el allocatable :: ihpb_all,jhpb_all !(maxss,maxstr_proc)
2096 ! allocate(icc(n_conf),iscore(n_conf)) !(maxconf)
2099 ! allocatable :: tempfac !(2,maxres)
2102 !el allocate(beta_h(maxT)) !(maxT)
2104 ! end subroutine alloc_clust_arrays
2105 !-----------------------------------------------------------------------------
2106 !-----------------------------------------------------------------------------