1 subroutine read_general_data(*)
4 include "DIMENSIONS.ZSCOPT"
5 include "DIMENSIONS.FREE"
6 include "COMMON.TORSION"
7 include "COMMON.INTERACT"
8 include "COMMON.IOUNITS"
11 include "COMMON.PROTFILES"
12 include "COMMON.CHAIN"
13 include "COMMON.NAMES"
14 include "COMMON.FFIELD"
15 include "COMMON.ENEPS"
16 include "COMMON.WEIGHTS"
18 include "COMMON.CONTROL"
19 include "COMMON.ENERGIES"
20 include "COMMON.SBRIDGE"
21 character*800 controlcard
22 integer i,j,k,ii,n_ene_found
23 integer ind,itype1,itype2,itypf,itypsc,itypp
30 call card_concat(controlcard,.true.)
31 call readi(controlcard,"N_ENE",n_ene,max_ene)
32 if (n_ene.gt.max_ene) then
33 write (iout,*) "Error: parameter out of range: N_ENE",n_ene,
37 call readi(controlcard,"NPARMSET",nparmset,1)
38 separate_parset = index(controlcard,"SEPARATE_PARSET").gt.0
39 call readi(controlcard,"IPARMPRINT",iparmprint,1)
40 write (iout,*) "PARMPRINT",iparmprint
41 if (nparmset.gt.max_parm) then
42 write (iout,*) "Error: parameter out of range: NPARMSET",
46 call readi(controlcard,"MAXIT",maxit,5000)
47 call reada(controlcard,"FIMIN",fimin,1.0d-3)
48 call readi(controlcard,"ENSEMBLES",ensembles,0)
49 hamil_rep=index(controlcard,"HAMIL_REP").gt.0
50 write (iout,*) "Number of energy parameter sets",nparmset
51 call multreadi(controlcard,"ISAMPL",isampl,nparmset,1)
52 write (iout,*) "MaxSlice",MaxSlice
53 call readi(controlcard,"NSLICE",nslice,1)
54 call reada(controlcard,"FRAC_CUTOFF",frac_cutoff,10.0d0)
55 write (iout,*) "RMSD cutoff for fraction calculations",frac_cutoff
57 if (nslice.gt.MaxSlice) then
58 write (iout,*) "Error: parameter out of range: NSLICE",nslice,
62 write (iout,*) "Frequency of storing conformations",
63 & (isampl(i),i=1,nparmset)
64 write (iout,*) "Maxit",maxit," Fimin",fimin
65 call readi(controlcard,"NQ",nQ,1)
67 write (iout,*) "Error: parameter out of range: NQ",nq,
72 if (index(controlcard,"CLASSIFY").gt.0) indpdb=1
73 call reada(controlcard,"DELTA",delta,1.0d-2)
74 call readi(controlcard,"EINICHECK",einicheck,2)
75 call reada(controlcard,"DELTRMS",deltrms,5.0d-2)
76 call reada(controlcard,"DELTRGY",deltrgy,5.0d-2)
77 call readi(controlcard,"RESCALE",rescale_mode,1)
78 check_conf=index(controlcard,"NO_CHECK_CONF").eq.0
79 call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
80 call readi(controlcard,'SYM',symetr,1)
81 write (iout,*) "DISTCHAINMAX",distchainmax
82 write (iout,*) "delta",delta
83 write (iout,*) "einicheck",einicheck
84 write (iout,*) "rescale_mode",rescale_mode
86 bxfile=index(controlcard,"BXFILE").gt.0
87 cxfile=index(controlcard,"CXFILE").gt.0
88 if (nslice .eq. 1 .and. .not.bxfile .and. .not.cxfile)
90 histfile=index(controlcard,"HISTFILE").gt.0
91 histout=index(controlcard,"HISTOUT").gt.0
92 entfile=index(controlcard,"ENTFILE").gt.0
93 zscfile=index(controlcard,"ZSCFILE").gt.0
94 with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
95 write (iout,*) "with_dihed_constr ",with_dihed_constr
96 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
97 write (iout,*) "with_theta_constr ",with_theta_constr
98 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
99 dyn_ss=index(controlcard,"DYN_SS").gt.0
102 c------------------------------------------------------------------------------
103 subroutine read_efree(*)
105 C Read molecular data
109 include 'DIMENSIONS.ZSCOPT'
110 include 'DIMENSIONS.COMPAR'
111 include 'DIMENSIONS.FREE'
112 include 'COMMON.IOUNITS'
113 include 'COMMON.TIME1'
114 include 'COMMON.SBRIDGE'
115 include 'COMMON.CONTROL'
116 include 'COMMON.CHAIN'
117 include 'COMMON.HEADER'
119 include 'COMMON.FREE'
120 character*320 controlcard,ucase
121 integer iparm,ib,i,j,npars
133 call card_concat(controlcard,.true.)
134 call readi(controlcard,'NT',nT_h(iparm),1)
135 write (iout,*) "iparm",iparm," nt",nT_h(iparm)
137 if (nT_h(iparm).gt.MaxT_h) then
138 write (iout,*) "Error: parameter out of range: NT",nT_h(iparm),
142 replica(iparm)=index(controlcard,"REPLICA").gt.0
143 umbrella(iparm)=index(controlcard,"UMBRELLA").gt.0
144 read_iset(iparm)=index(controlcard,"READ_ISET").gt.0
145 write (iout,*) "nQ",nQ," nT",nT_h(iparm)," replica ",
146 & replica(iparm)," umbrella ",umbrella(iparm),
147 & " read_iset",read_iset(iparm)
150 call card_concat(controlcard,.true.)
151 call readi(controlcard,'NR',nR(ib,iparm),1)
152 if (umbrella(iparm)) then
155 nRR(ib,iparm)=nR(ib,iparm)
157 if (nR(ib,iparm).gt.MaxR) then
158 write (iout,*) "Error: parameter out of range: NR",
162 call reada(controlcard,'TEMP',beta_h(ib,iparm),298.0d0)
163 beta_h(ib,iparm)=1.0d0/(beta_h(ib,iparm)*1.987D-3)
164 call multreada(controlcard,'FI',f(1,ib,iparm),nR(ib,iparm),
167 call card_concat(controlcard,.true.)
168 call multreada(controlcard,'KH',KH(1,i,ib,iparm),nQ,
170 call multreada(controlcard,'Q0',Q0(1,i,ib,iparm),nQ,
175 write (iout,*) "ib",ib," beta_h",
176 & 1.0d0/(0.001987*beta_h(ib,iparm))
177 write (iout,*) "nR",nR(ib,iparm)
178 write (iout,*) "fi",(f(i,ib,iparm),i=1,nR(ib,iparm))
180 write (iout,*) "i",i," Kh",(Kh(j,i,ib,iparm),j=1,nQ),
181 & "q0",(q0(j,i,ib,iparm),j=1,nQ)
193 nR(ib,iparm)=nR(ib,1)
194 if (umbrella(iparm)) then
197 nRR(ib,iparm)=nR(ib,1)
199 beta_h(ib,iparm)=beta_h(ib,1)
201 f(i,ib,iparm)=f(i,ib,1)
203 KH(j,i,ib,iparm)=KH(j,i,ib,1)
204 Q0(j,i,ib,iparm)=Q0(j,i,ib,1)
207 replica(iparm)=replica(1)
208 umbrella(iparm)=umbrella(1)
209 read_iset(iparm)=read_iset(1)
217 c-----------------------------------------------------------------------------
218 subroutine read_protein_data(*)
221 include "DIMENSIONS.ZSCOPT"
222 include "DIMENSIONS.FREE"
225 integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
228 include "COMMON.CHAIN"
229 include "COMMON.IOUNITS"
230 include "COMMON.PROT"
231 include "COMMON.PROTFILES"
232 include "COMMON.NAMES"
233 include "COMMON.FREE"
234 include "COMMON.OBCINKA"
236 character*16000 controlcard
237 integer i,ii,ib,iR,iparm,ilen,iroof,nthr,npars
247 C Read names of files with conformation data.
248 if (replica(iparm)) then
254 do ii=1,nRR(ib,iparm)
255 write (iout,*) "Parameter set",iparm," temperature",ib,
258 call card_concat(controlcard,.true.)
259 write (iout,*) controlcard(:ilen(controlcard))
260 call readi(controlcard,"NFILE_BIN",nfile_bin(ii,ib,iparm),0)
261 call readi(controlcard,"NFILE_ASC",nfile_asc(ii,ib,iparm),0)
262 call readi(controlcard,"NFILE_CX",nfile_cx(ii,ib,iparm),0)
263 call readi(controlcard,"REC_START",rec_start(ii,ib,iparm),1)
264 call readi(controlcard,"REC_END",rec_end(ii,ib,iparm),
265 & maxstr*isampl(iparm)+rec_start(ii,ib,iparm)-1)
266 call reada(controlcard,"TIME_START",
267 & time_start_collect(ii,ib,iparm),0.0d0)
268 call reada(controlcard,"TIME_END",time_end_collect(ii,ib,iparm),
270 write (iout,*) "rec_start",rec_start(ii,ib,iparm),
271 & " rec_end",rec_end(ii,ib,iparm)
272 write (iout,*) "time_start",time_start_collect(ii,ib,iparm),
273 & " time_end",time_end_collect(ii,ib,iparm)
275 if (replica(iparm)) then
276 call readi(controlcard,"TOTRAJ",totraj(ii,iparm),1)
277 write (iout,*) "Number of trajectories",totraj(ii,iparm)
280 if (nfile_bin(ii,ib,iparm).lt.2
281 & .and. nfile_asc(ii,ib,iparm).eq.0
282 & .and. nfile_cx(ii,ib,iparm).eq.0) then
283 write (iout,*) "Error - no action specified!"
286 if (nfile_bin(ii,ib,iparm).gt.0) then
287 call card_concat(controlcard,.false.)
288 call split_string(controlcard,protfiles(1,1,ii,ib,iparm),
289 & maxfile_prot,nfile_bin(ii,ib,iparm))
291 write(iout,*)"nfile_bin",nfile_bin(ii,ib,iparm)
292 write(iout,*) (protfiles(i,1,ii,ib,iparm),
293 & i=1,nfile_bin(ii,ib,iparm))
296 if (nfile_asc(ii,ib,iparm).gt.0) then
297 call card_concat(controlcard,.false.)
298 call split_string(controlcard,protfiles(1,2,ii,ib,iparm),
299 & maxfile_prot,nfile_asc(ii,ib,iparm))
301 write(iout,*) "nfile_asc(ii,ib,iparm)",nfile_asc(ii,ib,iparm)
302 write(iout,*) (protfiles(i,2,ii,ib,iparm),
303 & i=1,nfile_asc(ii,ib,iparm))
305 else if (nfile_cx(ii,ib,iparm).gt.0) then
306 call card_concat(controlcard,.false.)
307 call split_string(controlcard,protfiles(1,2,ii,ib,iparm),
308 & maxfile_prot,nfile_cx(ii,ib,iparm))
310 write(iout,*) "nfile_cx(ii,ib,iparm)",nfile_cx(ii,ib,iparm)
311 write(iout,*) (protfiles(i,2,ii,ib,iparm),
312 & i=1,nfile_cx(ii,ib,iparm))
323 c-------------------------------------------------------------------------------
324 subroutine opentmp(islice,iunit,bprotfile_temp)
327 include "DIMENSIONS.ZSCOPT"
328 include "DIMENSIONS.FREE"
331 integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
334 include "COMMON.IOUNITS"
335 include "COMMON.PROTFILES"
336 include "COMMON.PROT"
337 include "COMMON.FREE"
338 character*64 bprotfile_temp
339 character*3 liczba,liczba2
346 write (liczba1,'(bz,i2.2)') islice
347 write (liczba,'(bz,i3.3)') me
349 c write (iout,*) "separate_parset ",separate_parset,
351 if (separate_parset) then
352 write (liczba2,'(bz,i3.3)') myparm
353 bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
354 & prefix(:ilen(prefix))//liczba//"_"//liczba2//".xbin.tmp"//liczba1
355 open (iunit,file=bprotfile_temp,status="unknown",
356 & form="unformatted",access="direct",recl=lenrec)
358 bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
359 & prefix(:ilen(prefix))//liczba//".xbin.tmp"//liczba1
360 open (iunit,file=bprotfile_temp,status="unknown",
361 & form="unformatted",access="direct",recl=lenrec)
364 bprotfile_temp = scratchdir(:ilen(scratchdir))//
365 & "/"//prefix(:ilen(prefix))//".xbin.tmp"//liczba1
366 open (iunit,file=bprotfile_temp,status="unknown",
367 & form="unformatted",access="direct",recl=lenrec)
369 c write (iout,*) "OpenTmp iunit",iunit," bprotfile_temp",
374 c-------------------------------------------------------------------------------
375 subroutine read_database(*)
378 include "DIMENSIONS.ZSCOPT"
379 include "DIMENSIONS.FREE"
382 integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
385 include "COMMON.CHAIN"
386 include "COMMON.IOUNITS"
387 include "COMMON.PROTFILES"
388 include "COMMON.NAMES"
391 include "COMMON.ENEPS"
392 include "COMMON.PROT"
393 include "COMMON.INTERACT"
394 include "COMMON.FREE"
395 include "COMMON.SBRIDGE"
396 include "COMMON.OBCINKA"
397 real*4 csingle(3,maxres2)
398 character*64 nazwa,bprotfile_temp
401 integer i,j,ii,jj(maxslice),k,kk(maxslice),l,
402 & ll(maxslice),mm(maxslice),if
403 integer nrec,nlines,iscor,iunit,islice
404 double precision energ
407 double precision rmsdev,energia(0:max_ene),efree,eini,temp
408 double precision prop(maxQ)
409 integer ntot_all(maxslice,0:maxprocs-1), maxslice_buff
410 integer iparm,ib,iib,ir,nprop,nthr,npars
411 double precision etot,time
415 lenrec1=12*(nres+nct-nnt+1)+4*(2*nss+2)+24
416 lenrec2=12*(nres+nct-nnt+1)+4*(2*nss+2)+24+8*nQ
418 write (iout,*) "lenrec",lenrec," lenrec1",lenrec1,
428 write (iout,*) "nparmset",nparmset
436 if (replica(iparm)) then
443 do iR=1,nRR(ib,iparm)
445 write (iout,*) "iparm",iparm," ib",ib," iR",iR," nQ",nQ
451 IF (NFILE_BIN(iR,ib,iparm).GT.0) THEN
452 c Read conformations from binary DA files (one per batch) and write them to
453 c a binary DA scratchfile.
454 write (liczba,'(bz,i3.3)') me
455 do if=1,nfile_bin(iR,ib,iparm)
456 nazwa=protfiles(if,1,iR,ib,iparm)
457 & (:ilen(protfiles(if,1,iR,ib,iparm)))//".bx"
458 open (ientin,file=nazwa,status="old",form="unformatted",
459 & access="direct",recl=lenrec2,err=1111)
462 call opentmp(islice,ientout,bprotfile_temp)
463 call bxread(nazwa,ii,jj(islice),kk(islice),ll(islice),
464 & mm(islice),iR,ib,iparm)
471 IF (NFILE_ASC(iR,ib,iparm).GT.0) THEN
472 c Read conformations from multiple ASCII int files and write them to a binary
474 do if=1,nfile_asc(iR,ib,iparm)
475 nazwa=protfiles(if,2,iR,ib,iparm)
476 & (:ilen(protfiles(if,2,iR,ib,iparm)))//".x"
477 open(unit=ientin,file=nazwa,status='old',err=1111)
478 write(iout,*) "reading ",nazwa(:ilen(nazwa))
480 call xread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm)
483 IF (NFILE_CX(iR,ib,iparm).gt.0) THEN
484 c Read conformations from cx files and write them to a binary
486 do if=1,nfile_cx(iR,ib,iparm)
487 nazwa=protfiles(if,2,iR,ib,iparm)
488 & (:ilen(protfiles(if,2,iR,ib,iparm)))//".cx"
489 write(iout,*) "reading ",nazwa(:ilen(nazwa))
491 print *,"Calling cxread"
492 call cxread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm,
495 write (iout,*) "exit cxread"
501 stot(islice)=stot(islice)+jj(islice)
506 write (iout,*) "IPARM",iparm
509 if (nslice.eq.1) then
511 write (liczba,'(bz,i3.3)') me
512 bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
513 & prefix(:ilen(prefix))//liczba//".xbin.tmp"
515 bprotfile_temp = scratchdir(:ilen(scratchdir))//
516 & "/"//prefix(:ilen(prefix))//".xbin.tmp"
518 write(iout,*) mm(1)," conformations read",ll(1),
519 & " conformations written to ",
520 & bprotfile_temp(:ilen(bprotfile_temp))
523 write (liczba1,'(bz,i2.2)') islice
525 write (liczba,'(bz,i3.3)') me
526 bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
527 & prefix(:ilen(prefix))//liczba//".xbin.tmp"//liczba1
529 bprotfile_temp = scratchdir(:ilen(scratchdir))//
530 & "/"//prefix(:ilen(prefix))//".xbin.tmp"//liczba1
532 write(iout,*) mm(islice)," conformations read",ll(islice),
533 & " conformations written to ",
534 & bprotfile_temp(:ilen(bprotfile_temp))
539 c Check if everyone has the same number of conformations
541 c call MPI_ALLgather(MPI_IN_PLACE,stot(1),MPI_DATATYPE_NULL,
542 c & ntot_all(1,0),maxslice,MPI_INTEGER,MPI_Comm_World,IERROR)
544 maxslice_buff=maxslice
546 call MPI_Allgather(stot(1),maxslice_buff,MPI_INTEGER,
547 & ntot_all(1,0),maxslice,MPI_INTEGER,MPI_Comm_World,IERROR)
552 if (stot(islice).ne.ntot_all(islice,i)) then
553 write (iout,*) "Number of conformations at processor",i,
554 & " differs from that at processor",me,
555 & stot(islice),ntot_all(islice,i)," slice",islice
563 write (iout,*) "Numbers of conformations read by processors"
566 write (iout,'(8i10)') i,(ntot_all(islice,i),islice=1,nslice)
568 write (iout,*) "Calculation terminated."
573 ntot(islice)=stot(islice)
577 1111 write(iout,*) "Error opening coordinate file ",nazwa(:ilen(nazwa))
581 c------------------------------------------------------------------------------
582 subroutine card_concat(card,to_upper)
584 include 'DIMENSIONS.ZSCOPT'
585 include "COMMON.IOUNITS"
587 character*80 karta,ucase
591 read (inp,'(a)') karta
592 if (to_upper) karta=ucase(karta)
594 do while (karta(80:80).eq.'&')
595 card=card(:ilen(card)+1)//karta(:79)
596 read (inp,'(a)') karta
597 if (to_upper) karta=ucase(karta)
599 card=card(:ilen(card)+1)//karta
602 c------------------------------------------------------------------------------
603 subroutine readi(rekord,lancuch,wartosc,default)
605 character*(*) rekord,lancuch
606 integer wartosc,default
609 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
614 iread=iread+ilen(lancuch)+1
615 read (rekord(iread:),*) wartosc
618 c----------------------------------------------------------------------------
619 subroutine reada(rekord,lancuch,wartosc,default)
621 character*(*) rekord,lancuch
623 double precision wartosc,default
626 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
631 iread=iread+ilen(lancuch)+1
632 read (rekord(iread:),*) wartosc
635 c----------------------------------------------------------------------------
636 subroutine multreadi(rekord,lancuch,tablica,dim,default)
639 integer tablica(dim),default
640 character*(*) rekord,lancuch
647 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
648 if (iread.eq.0) return
649 iread=iread+ilen(lancuch)+1
650 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
653 c----------------------------------------------------------------------------
654 subroutine multreada(rekord,lancuch,tablica,dim,default)
657 double precision tablica(dim),default
658 character*(*) rekord,lancuch
665 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
666 if (iread.eq.0) return
667 iread=iread+ilen(lancuch)+1
668 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
671 c----------------------------------------------------------------------------
672 subroutine reads(rekord,lancuch,wartosc,default)
674 character*(*) rekord,lancuch,wartosc,default
676 integer ilen,lenlan,lenrec,iread,ireade
682 iread=index(rekord,lancuch(:lenlan)//"=")
683 c print *,"rekord",rekord," lancuch",lancuch
684 c print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
690 c print *,"iread",iread
691 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
692 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
694 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
696 c print *,"iread",iread
697 if (iread.gt.lenrec) then
702 c print *,"ireade",ireade
703 do while (ireade.lt.lenrec .and.
704 & .not.iblnk(rekord(ireade:ireade)))
707 wartosc=rekord(iread:ireade)
710 c----------------------------------------------------------------------------
711 subroutine multreads(rekord,lancuch,tablica,dim,default)
714 character*(*) rekord,lancuch,tablica(dim),default
716 integer ilen,lenlan,lenrec,iread,ireade
725 iread=index(rekord,lancuch(:lenlan)//"=")
726 c print *,"rekord",rekord," lancuch",lancuch
727 c print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
728 if (iread.eq.0) return
731 c print *,"iread",iread
732 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
733 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
735 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
737 c print *,"iread",iread
738 if (iread.gt.lenrec) return
740 c print *,"ireade",ireade
741 do while (ireade.lt.lenrec .and.
742 & .not.iblnk(rekord(ireade:ireade)))
745 tablica(i)=rekord(iread:ireade)
749 c----------------------------------------------------------------------------
750 subroutine split_string(rekord,tablica,dim,nsub)
752 integer dim,nsub,i,ii,ll,kk
753 character*(*) tablica(dim)
764 C Find the start of term name
766 do while (ii.le.ll .and. rekord(ii:ii).eq." ")
769 C Parse the name into TABLICA(i) until blank found
770 do while (ii.le.ll .and. rekord(ii:ii).ne." ")
772 tablica(i)(kk:kk)=rekord(ii:ii)
775 if (kk.gt.0) nsub=nsub+1
780 c--------------------------------------------------------------------------------
781 integer function iroof(n,m)
783 if (ii*m .lt. n) ii=ii+1