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 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)
55 if (nslice.gt.MaxSlice) then
56 write (iout,*) "Error: parameter out of range: NSLICE",nslice,
60 write (iout,*) "Frequency of storing conformations",
61 & (isampl(i),i=1,nparmset)
62 write (iout,*) "Maxit",maxit," Fimin",fimin
63 call readi(controlcard,"NQ",nQ,1)
65 write (iout,*) "Error: parameter out of range: NQ",nq,
70 if (index(controlcard,"CLASSIFY").gt.0) indpdb=1
71 call reada(controlcard,"DELTA",delta,1.0d-2)
72 call readi(controlcard,"EINICHECK",einicheck,2)
73 call reada(controlcard,"DELTRMS",deltrms,5.0d-2)
74 call reada(controlcard,"DELTRGY",deltrgy,5.0d-2)
75 call readi(controlcard,"RESCALE",rescale_mode,1)
76 check_conf=index(controlcard,"NO_CHECK_CONF").eq.0
77 call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
78 call reada(controlcard,'BOXX',boxxsize,100.0d0)
79 call reada(controlcard,'BOXY',boxysize,100.0d0)
80 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
81 c Cutoff range for interactions
82 call reada(controlcard,"R_CUT",r_cut,15.0d0)
83 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
84 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
85 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
86 if (lipthick.gt.0.0d0) then
87 bordliptop=(boxzsize+lipthick)/2.0
88 bordlipbot=bordliptop-lipthick
90 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
91 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
92 buflipbot=bordlipbot+lipbufthick
93 bufliptop=bordliptop-lipbufthick
94 if ((lipbufthick*2.0d0).gt.lipthick)
95 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
97 write(iout,*) "bordliptop=",bordliptop
98 write(iout,*) "bordlipbot=",bordlipbot
99 write(iout,*) "bufliptop=",bufliptop
100 write(iout,*) "buflipbot=",buflipbot
101 call readi(controlcard,'SYM',symetr,1)
102 write (iout,*) "DISTCHAINMAX",distchainmax
103 write (iout,*) "delta",delta
104 write (iout,*) "einicheck",einicheck
105 write (iout,*) "rescale_mode",rescale_mode
107 bxfile=index(controlcard,"BXFILE").gt.0
108 cxfile=index(controlcard,"CXFILE").gt.0
109 if (nslice .eq. 1 .and. .not.bxfile .and. .not.cxfile)
111 histfile=index(controlcard,"HISTFILE").gt.0
112 histout=index(controlcard,"HISTOUT").gt.0
113 entfile=index(controlcard,"ENTFILE").gt.0
114 zscfile=index(controlcard,"ZSCFILE").gt.0
115 with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
116 write (iout,*) "with_dihed_constr ",with_dihed_constr
117 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
120 c------------------------------------------------------------------------------
121 subroutine read_efree(*)
123 C Read molecular data
127 include 'DIMENSIONS.ZSCOPT'
128 include 'DIMENSIONS.COMPAR'
129 include 'DIMENSIONS.FREE'
130 include 'COMMON.IOUNITS'
131 include 'COMMON.TIME1'
132 include 'COMMON.SBRIDGE'
133 include 'COMMON.CONTROL'
134 include 'COMMON.CHAIN'
135 include 'COMMON.HEADER'
137 include 'COMMON.FREE'
138 character*320 controlcard,ucase
139 integer iparm,ib,i,j,npars
151 call card_concat(controlcard,.true.)
152 call readi(controlcard,'NT',nT_h(iparm),1)
153 write (iout,*) "iparm",iparm," nt",nT_h(iparm)
155 if (nT_h(iparm).gt.MaxT_h) then
156 write (iout,*) "Error: parameter out of range: NT",nT_h(iparm),
160 replica(iparm)=index(controlcard,"REPLICA").gt.0
161 umbrella(iparm)=index(controlcard,"UMBRELLA").gt.0
162 read_iset(iparm)=index(controlcard,"READ_ISET").gt.0
163 write (iout,*) "nQ",nQ," nT",nT_h(iparm)," replica ",
164 & replica(iparm)," umbrella ",umbrella(iparm),
165 & " read_iset",read_iset(iparm)
168 call card_concat(controlcard,.true.)
169 call readi(controlcard,'NR',nR(ib,iparm),1)
170 if (umbrella(iparm)) then
173 nRR(ib,iparm)=nR(ib,iparm)
175 if (nR(ib,iparm).gt.MaxR) then
176 write (iout,*) "Error: parameter out of range: NR",
180 call reada(controlcard,'TEMP',beta_h(ib,iparm),298.0d0)
181 beta_h(ib,iparm)=1.0d0/(beta_h(ib,iparm)*1.987D-3)
182 call multreada(controlcard,'FI',f(1,ib,iparm),nR(ib,iparm),
185 call card_concat(controlcard,.true.)
186 call multreada(controlcard,'KH',KH(1,i,ib,iparm),nQ,
188 call multreada(controlcard,'Q0',Q0(1,i,ib,iparm),nQ,
193 write (iout,*) "ib",ib," beta_h",
194 & 1.0d0/(0.001987*beta_h(ib,iparm))
195 write (iout,*) "nR",nR(ib,iparm)
196 write (iout,*) "fi",(f(i,ib,iparm),i=1,nR(ib,iparm))
198 write (iout,*) "i",i," Kh",(Kh(j,i,ib,iparm),j=1,nQ),
199 & "q0",(q0(j,i,ib,iparm),j=1,nQ)
211 nR(ib,iparm)=nR(ib,1)
212 if (umbrella(iparm)) then
215 nRR(ib,iparm)=nR(ib,1)
217 beta_h(ib,iparm)=beta_h(ib,1)
219 f(i,ib,iparm)=f(i,ib,1)
221 KH(j,i,ib,iparm)=KH(j,i,ib,1)
222 Q0(j,i,ib,iparm)=Q0(j,i,ib,1)
225 replica(iparm)=replica(1)
226 umbrella(iparm)=umbrella(1)
227 read_iset(iparm)=read_iset(1)
235 c-----------------------------------------------------------------------------
236 subroutine read_protein_data(*)
239 include "DIMENSIONS.ZSCOPT"
240 include "DIMENSIONS.FREE"
243 integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
246 include "COMMON.CHAIN"
247 include "COMMON.IOUNITS"
248 include "COMMON.PROT"
249 include "COMMON.PROTFILES"
250 include "COMMON.NAMES"
251 include "COMMON.FREE"
252 include "COMMON.OBCINKA"
254 character*16000 controlcard
255 integer i,ii,ib,iR,iparm,ilen,iroof,nthr,npars
265 C Read names of files with conformation data.
266 if (replica(iparm)) then
272 do ii=1,nRR(ib,iparm)
273 write (iout,*) "Parameter set",iparm," temperature",ib,
276 call card_concat(controlcard,.true.)
277 write (iout,*) controlcard(:ilen(controlcard))
278 call readi(controlcard,"NFILE_BIN",nfile_bin(ii,ib,iparm),0)
279 call readi(controlcard,"NFILE_ASC",nfile_asc(ii,ib,iparm),0)
280 call readi(controlcard,"NFILE_CX",nfile_cx(ii,ib,iparm),0)
281 call readi(controlcard,"REC_START",rec_start(ii,ib,iparm),1)
282 call readi(controlcard,"REC_END",rec_end(ii,ib,iparm),
283 & maxstr*isampl(iparm)+rec_start(ii,ib,iparm)-1)
284 call reada(controlcard,"TIME_START",
285 & time_start_collect(ii,ib,iparm),0.0d0)
286 call reada(controlcard,"TIME_END",time_end_collect(ii,ib,iparm),
288 write (iout,*) "rec_start",rec_start(ii,ib,iparm),
289 & " rec_end",rec_end(ii,ib,iparm)
290 write (iout,*) "time_start",time_start_collect(ii,ib,iparm),
291 & " time_end",time_end_collect(ii,ib,iparm)
293 if (replica(iparm)) then
294 call readi(controlcard,"TOTRAJ",totraj(ii,iparm),1)
295 write (iout,*) "Number of trajectories",totraj(ii,iparm)
298 if (nfile_bin(ii,ib,iparm).lt.2
299 & .and. nfile_asc(ii,ib,iparm).eq.0
300 & .and. nfile_cx(ii,ib,iparm).eq.0) then
301 write (iout,*) "Error - no action specified!"
304 if (nfile_bin(ii,ib,iparm).gt.0) then
305 call card_concat(controlcard,.false.)
306 call split_string(controlcard,protfiles(1,1,ii,ib,iparm),
307 & maxfile_prot,nfile_bin(ii,ib,iparm))
309 write(iout,*)"nfile_bin",nfile_bin(ii,ib,iparm)
310 write(iout,*) (protfiles(i,1,ii,ib,iparm),
311 & i=1,nfile_bin(ii,ib,iparm))
314 if (nfile_asc(ii,ib,iparm).gt.0) then
315 call card_concat(controlcard,.false.)
316 call split_string(controlcard,protfiles(1,2,ii,ib,iparm),
317 & maxfile_prot,nfile_asc(ii,ib,iparm))
319 write(iout,*) "nfile_asc(ii,ib,iparm)",nfile_asc(ii,ib,iparm)
320 write(iout,*) (protfiles(i,2,ii,ib,iparm),
321 & i=1,nfile_asc(ii,ib,iparm))
323 else if (nfile_cx(ii,ib,iparm).gt.0) then
324 call card_concat(controlcard,.false.)
325 call split_string(controlcard,protfiles(1,2,ii,ib,iparm),
326 & maxfile_prot,nfile_cx(ii,ib,iparm))
328 write(iout,*) "nfile_cx(ii,ib,iparm)",nfile_cx(ii,ib,iparm)
329 write(iout,*) (protfiles(i,2,ii,ib,iparm),
330 & i=1,nfile_cx(ii,ib,iparm))
341 c-------------------------------------------------------------------------------
342 subroutine opentmp(islice,iunit,bprotfile_temp)
345 include "DIMENSIONS.ZSCOPT"
346 include "DIMENSIONS.FREE"
349 integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
352 include "COMMON.IOUNITS"
353 include "COMMON.PROTFILES"
354 include "COMMON.PROT"
355 include "COMMON.FREE"
356 character*64 bprotfile_temp
357 character*3 liczba,liczba2
364 write (liczba1,'(bz,i2.2)') islice
365 write (liczba,'(bz,i3.3)') me
367 c write (iout,*) "separate_parset ",separate_parset,
369 if (separate_parset) then
370 write (liczba2,'(bz,i3.3)') myparm
371 bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
372 & prefix(:ilen(prefix))//liczba//"_"//liczba2//".xbin.tmp"//liczba1
373 open (iunit,file=bprotfile_temp,status="unknown",
374 & form="unformatted",access="direct",recl=lenrec)
376 bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
377 & prefix(:ilen(prefix))//liczba//".xbin.tmp"//liczba1
378 open (iunit,file=bprotfile_temp,status="unknown",
379 & form="unformatted",access="direct",recl=lenrec)
382 bprotfile_temp = scratchdir(:ilen(scratchdir))//
383 & "/"//prefix(:ilen(prefix))//".xbin.tmp"//liczba1
384 open (iunit,file=bprotfile_temp,status="unknown",
385 & form="unformatted",access="direct",recl=lenrec)
387 c write (iout,*) "OpenTmp iunit",iunit," bprotfile_temp",
392 c-------------------------------------------------------------------------------
393 subroutine read_database(*)
396 include "DIMENSIONS.ZSCOPT"
397 include "DIMENSIONS.FREE"
400 integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
403 include "COMMON.CHAIN"
404 include "COMMON.IOUNITS"
405 include "COMMON.PROTFILES"
406 include "COMMON.NAMES"
409 include "COMMON.ENEPS"
410 include "COMMON.PROT"
411 include "COMMON.INTERACT"
412 include "COMMON.FREE"
413 include "COMMON.SBRIDGE"
414 include "COMMON.OBCINKA"
415 real*4 csingle(3,maxres2)
416 character*64 nazwa,bprotfile_temp
419 integer i,j,ii,jj(maxslice),k,kk(maxslice),l,
420 & ll(maxslice),mm(maxslice),if
421 integer nrec,nlines,iscor,iunit,islice
422 double precision energ
425 double precision rmsdev,energia(0:max_ene),efree,eini,temp
426 double precision prop(maxQ)
427 integer ntot_all(maxslice,0:maxprocs-1)
428 integer iparm,ib,iib,ir,nprop,nthr,npars
429 double precision etot,time
433 lenrec1=12*(nres+nct-nnt+1)+4*(2*nss+2)+24
434 lenrec2=12*(nres+nct-nnt+1)+4*(2*nss+2)+24+8*nQ
436 write (iout,*) "lenrec",lenrec," lenrec1",lenrec1,
446 write (iout,*) "nparmset",nparmset
454 if (replica(iparm)) then
461 do iR=1,nRR(ib,iparm)
463 write (iout,*) "iparm",iparm," ib",ib," iR",iR," nQ",nQ
469 IF (NFILE_BIN(iR,ib,iparm).GT.0) THEN
470 c Read conformations from binary DA files (one per batch) and write them to
471 c a binary DA scratchfile.
472 write (liczba,'(bz,i3.3)') me
473 do if=1,nfile_bin(iR,ib,iparm)
474 nazwa=protfiles(if,1,iR,ib,iparm)
475 & (:ilen(protfiles(if,1,iR,ib,iparm)))//".bx"
476 open (ientin,file=nazwa,status="old",form="unformatted",
477 & access="direct",recl=lenrec2,err=1111)
480 call opentmp(islice,ientout,bprotfile_temp)
481 call bxread(nazwa,ii,jj(islice),kk(islice),ll(islice),
482 & mm(islice),iR,ib,iparm)
489 IF (NFILE_ASC(iR,ib,iparm).GT.0) THEN
490 c Read conformations from multiple ASCII int files and write them to a binary
492 do if=1,nfile_asc(iR,ib,iparm)
493 nazwa=protfiles(if,2,iR,ib,iparm)
494 & (:ilen(protfiles(if,2,iR,ib,iparm)))//".x"
495 open(unit=ientin,file=nazwa,status='old',err=1111)
496 write(iout,*) "reading ",nazwa(:ilen(nazwa))
498 call xread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm)
501 IF (NFILE_CX(iR,ib,iparm).gt.0) THEN
502 c Read conformations from cx files and write them to a binary
504 do if=1,nfile_cx(iR,ib,iparm)
505 nazwa=protfiles(if,2,iR,ib,iparm)
506 & (:ilen(protfiles(if,2,iR,ib,iparm)))//".cx"
507 write(iout,*) "reading ",nazwa(:ilen(nazwa))
509 print *,"Calling cxread"
510 call cxread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm,
513 write (iout,*) "exit cxread"
519 stot(islice)=stot(islice)+jj(islice)
524 write (iout,*) "IPARM",iparm
527 if (nslice.eq.1) then
529 write (liczba,'(bz,i3.3)') me
530 bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
531 & prefix(:ilen(prefix))//liczba//".xbin.tmp"
533 bprotfile_temp = scratchdir(:ilen(scratchdir))//
534 & "/"//prefix(:ilen(prefix))//".xbin.tmp"
536 write(iout,*) mm(1)," conformations read",ll(1),
537 & " conformations written to ",
538 & bprotfile_temp(:ilen(bprotfile_temp))
541 write (liczba1,'(bz,i2.2)') islice
543 write (liczba,'(bz,i3.3)') me
544 bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
545 & prefix(:ilen(prefix))//liczba//".xbin.tmp"//liczba1
547 bprotfile_temp = scratchdir(:ilen(scratchdir))//
548 & "/"//prefix(:ilen(prefix))//".xbin.tmp"//liczba1
550 write(iout,*) mm(islice)," conformations read",ll(islice),
551 & " conformations written to ",
552 & bprotfile_temp(:ilen(bprotfile_temp))
557 c Check if everyone has the same number of conformations
558 call MPI_Allgather(stot(1),maxslice,MPI_INTEGER,
559 & ntot_all(1,0),maxslice,MPI_INTEGER,MPI_Comm_World,IERROR)
564 if (stot(islice).ne.ntot_all(islice,i)) then
565 write (iout,*) "Number of conformations at processor",i,
566 & " differs from that at processor",me,
567 & stot(islice),ntot_all(islice,i)," slice",islice
575 write (iout,*) "Numbers of conformations read by processors"
578 write (iout,'(8i10)') i,(ntot_all(islice,i),islice=1,nslice)
580 write (iout,*) "Calculation terminated."
585 ntot(islice)=stot(islice)
589 1111 write(iout,*) "Error opening coordinate file ",nazwa(:ilen(nazwa))
593 c------------------------------------------------------------------------------
594 subroutine card_concat(card,to_upper)
596 include 'DIMENSIONS.ZSCOPT'
597 include "COMMON.IOUNITS"
599 character*80 karta,ucase
603 read (inp,'(a)') karta
604 if (to_upper) karta=ucase(karta)
606 do while (karta(80:80).eq.'&')
607 card=card(:ilen(card)+1)//karta(:79)
608 read (inp,'(a)') karta
609 if (to_upper) karta=ucase(karta)
611 card=card(:ilen(card)+1)//karta
614 c------------------------------------------------------------------------------
615 subroutine readi(rekord,lancuch,wartosc,default)
617 character*(*) rekord,lancuch
618 integer wartosc,default
621 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
626 iread=iread+ilen(lancuch)+1
627 read (rekord(iread:),*) wartosc
630 c----------------------------------------------------------------------------
631 subroutine reada(rekord,lancuch,wartosc,default)
633 character*(*) rekord,lancuch
635 double precision wartosc,default
638 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
643 iread=iread+ilen(lancuch)+1
644 read (rekord(iread:),*) wartosc
647 c----------------------------------------------------------------------------
648 subroutine multreadi(rekord,lancuch,tablica,dim,default)
651 integer tablica(dim),default
652 character*(*) rekord,lancuch
659 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
660 if (iread.eq.0) return
661 iread=iread+ilen(lancuch)+1
662 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
665 c----------------------------------------------------------------------------
666 subroutine multreada(rekord,lancuch,tablica,dim,default)
669 double precision tablica(dim),default
670 character*(*) rekord,lancuch
677 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
678 if (iread.eq.0) return
679 iread=iread+ilen(lancuch)+1
680 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
683 c----------------------------------------------------------------------------
684 subroutine reads(rekord,lancuch,wartosc,default)
686 character*(*) rekord,lancuch,wartosc,default
688 integer ilen,lenlan,lenrec,iread,ireade
694 iread=index(rekord,lancuch(:lenlan)//"=")
695 c print *,"rekord",rekord," lancuch",lancuch
696 c print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
702 c print *,"iread",iread
703 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
704 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
706 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
708 c print *,"iread",iread
709 if (iread.gt.lenrec) then
714 c print *,"ireade",ireade
715 do while (ireade.lt.lenrec .and.
716 & .not.iblnk(rekord(ireade:ireade)))
719 wartosc=rekord(iread:ireade)
722 c----------------------------------------------------------------------------
723 subroutine multreads(rekord,lancuch,tablica,dim,default)
726 character*(*) rekord,lancuch,tablica(dim),default
728 integer ilen,lenlan,lenrec,iread,ireade
737 iread=index(rekord,lancuch(:lenlan)//"=")
738 c print *,"rekord",rekord," lancuch",lancuch
739 c print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
740 if (iread.eq.0) return
743 c print *,"iread",iread
744 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
745 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
747 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
749 c print *,"iread",iread
750 if (iread.gt.lenrec) return
752 c print *,"ireade",ireade
753 do while (ireade.lt.lenrec .and.
754 & .not.iblnk(rekord(ireade:ireade)))
757 tablica(i)=rekord(iread:ireade)
761 c----------------------------------------------------------------------------
762 subroutine split_string(rekord,tablica,dim,nsub)
764 integer dim,nsub,i,ii,ll,kk
765 character*(*) tablica(dim)
776 C Find the start of term name
778 do while (ii.le.ll .and. rekord(ii:ii).eq." ")
781 C Parse the name into TABLICA(i) until blank found
782 do while (ii.le.ll .and. rekord(ii:ii).ne." ")
784 tablica(i)(kk:kk)=rekord(ii:ii)
787 if (kk.gt.0) nsub=nsub+1
792 c--------------------------------------------------------------------------------
793 integer function iroof(n,m)
795 if (ii*m .lt. n) ii=ii+1