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.SPLITELE"
21 include "COMMON.SBRIDGE"
22 character*800 controlcard
23 integer i,j,k,ii,n_ene_found
24 integer ind,itype1,itype2,itypf,itypsc,itypp
31 call card_concat(controlcard,.true.)
32 call readi(controlcard,"N_ENE",n_ene,max_ene)
33 if (n_ene.gt.max_ene) then
34 write (iout,*) "Error: parameter out of range: N_ENE",n_ene,
38 call readi(controlcard,"NPARMSET",nparmset,1)
39 separate_parset = index(controlcard,"SEPARATE_PARSET").gt.0
40 call readi(controlcard,"IPARMPRINT",iparmprint,1)
41 write (iout,*) "PARMPRINT",iparmprint
42 if (nparmset.gt.max_parm) then
43 write (iout,*) "Error: parameter out of range: NPARMSET",
47 call readi(controlcard,"MAXIT",maxit,5000)
48 call reada(controlcard,"FIMIN",fimin,1.0d-3)
49 call readi(controlcard,"ENSEMBLES",ensembles,0)
50 hamil_rep=index(controlcard,"HAMIL_REP").gt.0
51 write (iout,*) "Number of energy parameter sets",nparmset
52 call multreadi(controlcard,"ISAMPL",isampl,nparmset,1)
53 write (iout,*) "MaxSlice",MaxSlice
54 call readi(controlcard,"NSLICE",nslice,1)
56 if (nslice.gt.MaxSlice) then
57 write (iout,*) "Error: parameter out of range: NSLICE",nslice,
61 write (iout,*) "Frequency of storing conformations",
62 & (isampl(i),i=1,nparmset)
63 write (iout,*) "Maxit",maxit," Fimin",fimin
64 call readi(controlcard,"NQ",nQ,1)
66 write (iout,*) "Error: parameter out of range: NQ",nq,
71 if (index(controlcard,"CLASSIFY").gt.0) indpdb=1
72 call reada(controlcard,"DELTA",delta,1.0d-2)
73 call readi(controlcard,"EINICHECK",einicheck,2)
74 call reada(controlcard,"DELTRMS",deltrms,5.0d-2)
75 call reada(controlcard,"DELTRGY",deltrgy,5.0d-2)
76 call readi(controlcard,"RESCALE",rescale_mode,1)
77 check_conf=index(controlcard,"NO_CHECK_CONF").eq.0
78 call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
79 call reada(controlcard,'BOXX',boxxsize,100.0d0)
80 call reada(controlcard,'BOXY',boxysize,100.0d0)
81 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
82 c Cutoff range for interactions
83 call reada(controlcard,"R_CUT",r_cut,15.0d0)
84 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
85 call readi(controlcard,'SYM',symetr,1)
86 write (iout,*) "DISTCHAINMAX",distchainmax
87 write (iout,*) "delta",delta
88 write (iout,*) "einicheck",einicheck
89 write (iout,*) "rescale_mode",rescale_mode
91 bxfile=index(controlcard,"BXFILE").gt.0
92 cxfile=index(controlcard,"CXFILE").gt.0
93 if (nslice .eq. 1 .and. .not.bxfile .and. .not.cxfile)
95 histfile=index(controlcard,"HISTFILE").gt.0
96 histout=index(controlcard,"HISTOUT").gt.0
97 entfile=index(controlcard,"ENTFILE").gt.0
98 zscfile=index(controlcard,"ZSCFILE").gt.0
99 with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
100 write (iout,*) "with_dihed_constr ",with_dihed_constr
101 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
102 write (iout,*) "with_theta_constr ",with_theta_constr
103 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
104 dyn_ss=index(controlcard,"DYN_SS").gt.0
107 c------------------------------------------------------------------------------
108 subroutine read_efree(*)
110 C Read molecular data
114 include 'DIMENSIONS.ZSCOPT'
115 include 'DIMENSIONS.COMPAR'
116 include 'DIMENSIONS.FREE'
117 include 'COMMON.IOUNITS'
118 include 'COMMON.TIME1'
119 include 'COMMON.SBRIDGE'
120 include 'COMMON.CONTROL'
121 include 'COMMON.CHAIN'
122 include 'COMMON.HEADER'
124 include 'COMMON.FREE'
125 character*320 controlcard,ucase
126 integer iparm,ib,i,j,npars
138 call card_concat(controlcard,.true.)
139 call readi(controlcard,'NT',nT_h(iparm),1)
140 write (iout,*) "iparm",iparm," nt",nT_h(iparm)
142 if (nT_h(iparm).gt.MaxT_h) then
143 write (iout,*) "Error: parameter out of range: NT",nT_h(iparm),
147 replica(iparm)=index(controlcard,"REPLICA").gt.0
148 umbrella(iparm)=index(controlcard,"UMBRELLA").gt.0
149 read_iset(iparm)=index(controlcard,"READ_ISET").gt.0
150 write (iout,*) "nQ",nQ," nT",nT_h(iparm)," replica ",
151 & replica(iparm)," umbrella ",umbrella(iparm),
152 & " read_iset",read_iset(iparm)
155 call card_concat(controlcard,.true.)
156 call readi(controlcard,'NR',nR(ib,iparm),1)
157 if (umbrella(iparm)) then
160 nRR(ib,iparm)=nR(ib,iparm)
162 if (nR(ib,iparm).gt.MaxR) then
163 write (iout,*) "Error: parameter out of range: NR",
167 call reada(controlcard,'TEMP',beta_h(ib,iparm),298.0d0)
168 beta_h(ib,iparm)=1.0d0/(beta_h(ib,iparm)*1.987D-3)
169 call multreada(controlcard,'FI',f(1,ib,iparm),nR(ib,iparm),
172 call card_concat(controlcard,.true.)
173 call multreada(controlcard,'KH',KH(1,i,ib,iparm),nQ,
175 call multreada(controlcard,'Q0',Q0(1,i,ib,iparm),nQ,
180 write (iout,*) "ib",ib," beta_h",
181 & 1.0d0/(0.001987*beta_h(ib,iparm))
182 write (iout,*) "nR",nR(ib,iparm)
183 write (iout,*) "fi",(f(i,ib,iparm),i=1,nR(ib,iparm))
185 write (iout,*) "i",i," Kh",(Kh(j,i,ib,iparm),j=1,nQ),
186 & "q0",(q0(j,i,ib,iparm),j=1,nQ)
198 nR(ib,iparm)=nR(ib,1)
199 if (umbrella(iparm)) then
202 nRR(ib,iparm)=nR(ib,1)
204 beta_h(ib,iparm)=beta_h(ib,1)
206 f(i,ib,iparm)=f(i,ib,1)
208 KH(j,i,ib,iparm)=KH(j,i,ib,1)
209 Q0(j,i,ib,iparm)=Q0(j,i,ib,1)
212 replica(iparm)=replica(1)
213 umbrella(iparm)=umbrella(1)
214 read_iset(iparm)=read_iset(1)
222 c-----------------------------------------------------------------------------
223 subroutine read_protein_data(*)
226 include "DIMENSIONS.ZSCOPT"
227 include "DIMENSIONS.FREE"
230 integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
233 include "COMMON.CHAIN"
234 include "COMMON.IOUNITS"
235 include "COMMON.PROT"
236 include "COMMON.PROTFILES"
237 include "COMMON.NAMES"
238 include "COMMON.FREE"
239 include "COMMON.OBCINKA"
241 character*16000 controlcard
242 integer i,ii,ib,iR,iparm,ilen,iroof,nthr,npars
252 C Read names of files with conformation data.
253 if (replica(iparm)) then
259 do ii=1,nRR(ib,iparm)
260 write (iout,*) "Parameter set",iparm," temperature",ib,
263 call card_concat(controlcard,.true.)
264 write (iout,*) controlcard(:ilen(controlcard))
265 call readi(controlcard,"NFILE_BIN",nfile_bin(ii,ib,iparm),0)
266 call readi(controlcard,"NFILE_ASC",nfile_asc(ii,ib,iparm),0)
267 call readi(controlcard,"NFILE_CX",nfile_cx(ii,ib,iparm),0)
268 call readi(controlcard,"REC_START",rec_start(ii,ib,iparm),1)
269 call readi(controlcard,"REC_END",rec_end(ii,ib,iparm),
270 & maxstr*isampl(iparm)+rec_start(ii,ib,iparm)-1)
271 call reada(controlcard,"TIME_START",
272 & time_start_collect(ii,ib,iparm),0.0d0)
273 call reada(controlcard,"TIME_END",time_end_collect(ii,ib,iparm),
275 write (iout,*) "rec_start",rec_start(ii,ib,iparm),
276 & " rec_end",rec_end(ii,ib,iparm)
277 write (iout,*) "time_start",time_start_collect(ii,ib,iparm),
278 & " time_end",time_end_collect(ii,ib,iparm)
280 if (replica(iparm)) then
281 call readi(controlcard,"TOTRAJ",totraj(ii,iparm),1)
282 write (iout,*) "Number of trajectories",totraj(ii,iparm)
285 if (nfile_bin(ii,ib,iparm).lt.2
286 & .and. nfile_asc(ii,ib,iparm).eq.0
287 & .and. nfile_cx(ii,ib,iparm).eq.0) then
288 write (iout,*) "Error - no action specified!"
291 if (nfile_bin(ii,ib,iparm).gt.0) then
292 call card_concat(controlcard,.false.)
293 call split_string(controlcard,protfiles(1,1,ii,ib,iparm),
294 & maxfile_prot,nfile_bin(ii,ib,iparm))
296 write(iout,*)"nfile_bin",nfile_bin(ii,ib,iparm)
297 write(iout,*) (protfiles(i,1,ii,ib,iparm),
298 & i=1,nfile_bin(ii,ib,iparm))
301 if (nfile_asc(ii,ib,iparm).gt.0) then
302 call card_concat(controlcard,.false.)
303 call split_string(controlcard,protfiles(1,2,ii,ib,iparm),
304 & maxfile_prot,nfile_asc(ii,ib,iparm))
306 write(iout,*) "nfile_asc(ii,ib,iparm)",nfile_asc(ii,ib,iparm)
307 write(iout,*) (protfiles(i,2,ii,ib,iparm),
308 & i=1,nfile_asc(ii,ib,iparm))
310 else if (nfile_cx(ii,ib,iparm).gt.0) then
311 call card_concat(controlcard,.false.)
312 call split_string(controlcard,protfiles(1,2,ii,ib,iparm),
313 & maxfile_prot,nfile_cx(ii,ib,iparm))
315 write(iout,*) "nfile_cx(ii,ib,iparm)",nfile_cx(ii,ib,iparm)
316 write(iout,*) (protfiles(i,2,ii,ib,iparm),
317 & i=1,nfile_cx(ii,ib,iparm))
328 c-------------------------------------------------------------------------------
329 subroutine opentmp(islice,iunit,bprotfile_temp)
332 include "DIMENSIONS.ZSCOPT"
333 include "DIMENSIONS.FREE"
336 integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
339 include "COMMON.IOUNITS"
340 include "COMMON.PROTFILES"
341 include "COMMON.PROT"
342 include "COMMON.FREE"
343 character*64 bprotfile_temp
344 character*3 liczba,liczba2
351 write (liczba1,'(bz,i2.2)') islice
352 write (liczba,'(bz,i3.3)') me
354 c write (iout,*) "separate_parset ",separate_parset,
356 if (separate_parset) then
357 write (liczba2,'(bz,i3.3)') myparm
358 bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
359 & prefix(:ilen(prefix))//liczba//"_"//liczba2//".xbin.tmp"//liczba1
360 open (iunit,file=bprotfile_temp,status="unknown",
361 & form="unformatted",access="direct",recl=lenrec)
363 bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
364 & prefix(:ilen(prefix))//liczba//".xbin.tmp"//liczba1
365 open (iunit,file=bprotfile_temp,status="unknown",
366 & form="unformatted",access="direct",recl=lenrec)
369 bprotfile_temp = scratchdir(:ilen(scratchdir))//
370 & "/"//prefix(:ilen(prefix))//".xbin.tmp"//liczba1
371 open (iunit,file=bprotfile_temp,status="unknown",
372 & form="unformatted",access="direct",recl=lenrec)
374 c write (iout,*) "OpenTmp iunit",iunit," bprotfile_temp",
379 c-------------------------------------------------------------------------------
380 subroutine read_database(*)
383 include "DIMENSIONS.ZSCOPT"
384 include "DIMENSIONS.FREE"
387 integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
390 include "COMMON.CHAIN"
391 include "COMMON.IOUNITS"
392 include "COMMON.PROTFILES"
393 include "COMMON.NAMES"
396 include "COMMON.ENEPS"
397 include "COMMON.PROT"
398 include "COMMON.INTERACT"
399 include "COMMON.FREE"
400 include "COMMON.SBRIDGE"
401 include "COMMON.OBCINKA"
402 real*4 csingle(3,maxres2)
403 character*64 nazwa,bprotfile_temp
406 integer i,j,ii,jj(maxslice),k,kk(maxslice),l,
407 & ll(maxslice),mm(maxslice),if
408 integer nrec,nlines,iscor,iunit,islice
409 double precision energ
412 double precision rmsdev,energia(0:max_ene),efree,eini,temp
413 double precision prop(maxQ)
414 integer ntot_all(maxslice,0:maxprocs-1), maxslice_buff
415 integer iparm,ib,iib,ir,nprop,nthr,npars
416 double precision etot,time
420 lenrec1=12*(nres+nct-nnt+1)+4*(2*nss+2)+24
421 lenrec2=12*(nres+nct-nnt+1)+4*(2*nss+2)+24+8*nQ
423 write (iout,*) "lenrec",lenrec," lenrec1",lenrec1,
433 write (iout,*) "nparmset",nparmset
441 if (replica(iparm)) then
448 do iR=1,nRR(ib,iparm)
450 write (iout,*) "iparm",iparm," ib",ib," iR",iR," nQ",nQ
456 IF (NFILE_BIN(iR,ib,iparm).GT.0) THEN
457 c Read conformations from binary DA files (one per batch) and write them to
458 c a binary DA scratchfile.
459 write (liczba,'(bz,i3.3)') me
460 do if=1,nfile_bin(iR,ib,iparm)
461 nazwa=protfiles(if,1,iR,ib,iparm)
462 & (:ilen(protfiles(if,1,iR,ib,iparm)))//".bx"
463 open (ientin,file=nazwa,status="old",form="unformatted",
464 & access="direct",recl=lenrec2,err=1111)
467 call opentmp(islice,ientout,bprotfile_temp)
468 call bxread(nazwa,ii,jj(islice),kk(islice),ll(islice),
469 & mm(islice),iR,ib,iparm)
476 IF (NFILE_ASC(iR,ib,iparm).GT.0) THEN
477 c Read conformations from multiple ASCII int files and write them to a binary
479 do if=1,nfile_asc(iR,ib,iparm)
480 nazwa=protfiles(if,2,iR,ib,iparm)
481 & (:ilen(protfiles(if,2,iR,ib,iparm)))//".x"
482 open(unit=ientin,file=nazwa,status='old',err=1111)
483 write(iout,*) "reading ",nazwa(:ilen(nazwa))
485 call xread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm)
488 IF (NFILE_CX(iR,ib,iparm).gt.0) THEN
489 c Read conformations from cx files and write them to a binary
491 do if=1,nfile_cx(iR,ib,iparm)
492 nazwa=protfiles(if,2,iR,ib,iparm)
493 & (:ilen(protfiles(if,2,iR,ib,iparm)))//".cx"
494 write(iout,*) "reading ",nazwa(:ilen(nazwa))
496 print *,"Calling cxread"
497 call cxread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm,
500 write (iout,*) "exit cxread"
506 stot(islice)=stot(islice)+jj(islice)
511 write (iout,*) "IPARM",iparm
514 if (nslice.eq.1) then
516 write (liczba,'(bz,i3.3)') me
517 bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
518 & prefix(:ilen(prefix))//liczba//".xbin.tmp"
520 bprotfile_temp = scratchdir(:ilen(scratchdir))//
521 & "/"//prefix(:ilen(prefix))//".xbin.tmp"
523 write(iout,*) mm(1)," conformations read",ll(1),
524 & " conformations written to ",
525 & bprotfile_temp(:ilen(bprotfile_temp))
528 write (liczba1,'(bz,i2.2)') islice
530 write (liczba,'(bz,i3.3)') me
531 bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
532 & prefix(:ilen(prefix))//liczba//".xbin.tmp"//liczba1
534 bprotfile_temp = scratchdir(:ilen(scratchdir))//
535 & "/"//prefix(:ilen(prefix))//".xbin.tmp"//liczba1
537 write(iout,*) mm(islice)," conformations read",ll(islice),
538 & " conformations written to ",
539 & bprotfile_temp(:ilen(bprotfile_temp))
544 c Check if everyone has the same number of conformations
546 c call MPI_ALLgather(MPI_IN_PLACE,stot(1),MPI_DATATYPE_NULL,
547 c & ntot_all(1,0),maxslice,MPI_INTEGER,MPI_Comm_World,IERROR)
549 maxslice_buff=maxslice
551 call MPI_Allgather(stot(1),maxslice_buff,MPI_INTEGER,
552 & ntot_all(1,0),maxslice,MPI_INTEGER,MPI_Comm_World,IERROR)
557 if (stot(islice).ne.ntot_all(islice,i)) then
558 write (iout,*) "Number of conformations at processor",i,
559 & " differs from that at processor",me,
560 & stot(islice),ntot_all(islice,i)," slice",islice
568 write (iout,*) "Numbers of conformations read by processors"
571 write (iout,'(8i10)') i,(ntot_all(islice,i),islice=1,nslice)
573 write (iout,*) "Calculation terminated."
578 ntot(islice)=stot(islice)
582 1111 write(iout,*) "Error opening coordinate file ",nazwa(:ilen(nazwa))
586 c------------------------------------------------------------------------------
587 subroutine card_concat(card,to_upper)
589 include 'DIMENSIONS.ZSCOPT'
590 include "COMMON.IOUNITS"
592 character*80 karta,ucase
596 read (inp,'(a)') karta
597 if (to_upper) karta=ucase(karta)
599 do while (karta(80:80).eq.'&')
600 card=card(:ilen(card)+1)//karta(:79)
601 read (inp,'(a)') karta
602 if (to_upper) karta=ucase(karta)
604 card=card(:ilen(card)+1)//karta
607 c------------------------------------------------------------------------------
608 subroutine readi(rekord,lancuch,wartosc,default)
610 character*(*) rekord,lancuch
611 integer wartosc,default
614 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
619 iread=iread+ilen(lancuch)+1
620 read (rekord(iread:),*) wartosc
623 c----------------------------------------------------------------------------
624 subroutine reada(rekord,lancuch,wartosc,default)
626 character*(*) rekord,lancuch
628 double precision wartosc,default
631 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
636 iread=iread+ilen(lancuch)+1
637 read (rekord(iread:),*) wartosc
640 c----------------------------------------------------------------------------
641 subroutine multreadi(rekord,lancuch,tablica,dim,default)
644 integer tablica(dim),default
645 character*(*) rekord,lancuch
652 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
653 if (iread.eq.0) return
654 iread=iread+ilen(lancuch)+1
655 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
658 c----------------------------------------------------------------------------
659 subroutine multreada(rekord,lancuch,tablica,dim,default)
662 double precision tablica(dim),default
663 character*(*) rekord,lancuch
670 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
671 if (iread.eq.0) return
672 iread=iread+ilen(lancuch)+1
673 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
676 c----------------------------------------------------------------------------
677 subroutine reads(rekord,lancuch,wartosc,default)
679 character*(*) rekord,lancuch,wartosc,default
681 integer ilen,lenlan,lenrec,iread,ireade
687 iread=index(rekord,lancuch(:lenlan)//"=")
688 c print *,"rekord",rekord," lancuch",lancuch
689 c print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
695 c print *,"iread",iread
696 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
697 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
699 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
701 c print *,"iread",iread
702 if (iread.gt.lenrec) then
707 c print *,"ireade",ireade
708 do while (ireade.lt.lenrec .and.
709 & .not.iblnk(rekord(ireade:ireade)))
712 wartosc=rekord(iread:ireade)
715 c----------------------------------------------------------------------------
716 subroutine multreads(rekord,lancuch,tablica,dim,default)
719 character*(*) rekord,lancuch,tablica(dim),default
721 integer ilen,lenlan,lenrec,iread,ireade
730 iread=index(rekord,lancuch(:lenlan)//"=")
731 c print *,"rekord",rekord," lancuch",lancuch
732 c print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
733 if (iread.eq.0) return
736 c print *,"iread",iread
737 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
738 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
740 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
742 c print *,"iread",iread
743 if (iread.gt.lenrec) return
745 c print *,"ireade",ireade
746 do while (ireade.lt.lenrec .and.
747 & .not.iblnk(rekord(ireade:ireade)))
750 tablica(i)=rekord(iread:ireade)
754 c----------------------------------------------------------------------------
755 subroutine split_string(rekord,tablica,dim,nsub)
757 integer dim,nsub,i,ii,ll,kk
758 character*(*) tablica(dim)
769 C Find the start of term name
771 do while (ii.le.ll .and. rekord(ii:ii).eq." ")
774 C Parse the name into TABLICA(i) until blank found
775 do while (ii.le.ll .and. rekord(ii:ii).ne." ")
777 tablica(i)(kk:kk)=rekord(ii:ii)
780 if (kk.gt.0) nsub=nsub+1
785 c--------------------------------------------------------------------------------
786 integer function iroof(n,m)
788 if (ii*m .lt. n) ii=ii+1