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)
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 readi(controlcard,"NGRIDT",NGridT,400)
75 call reada(controlcard,"STARTGRIDT",StartGridT,200.0d0)
76 call reada(controlcard,"DELTA_T",Delta_T,1.0d0)
77 call reada(controlcard,"DELTRGY",deltrgy,5.0d-2)
78 call readi(controlcard,"RESCALE",rescale_mode,1)
79 check_conf=index(controlcard,"NO_CHECK_CONF").eq.0
80 write (iout,*) "delta",delta
81 write (iout,*) "einicheck",einicheck
82 write (iout,*) "rescale_mode",rescale_mode
84 bxfile=index(controlcard,"BXFILE").gt.0
85 cxfile=index(controlcard,"CXFILE").gt.0
86 if (nslice .eq. 1 .and. .not.bxfile .and. .not.cxfile)
88 histfile=index(controlcard,"HISTFILE").gt.0
89 histout=index(controlcard,"HISTOUT").gt.0
90 entfile=index(controlcard,"ENTFILE").gt.0
91 zscfile=index(controlcard,"ZSCFILE").gt.0
92 with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
93 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
94 write (iout,*) "with_dihed_constr ",with_dihed_constr,
95 & " CONSTR_DIST",constr_dist
96 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
97 write (iout,*) "with_homology_constr ",with_dihed_constr,
98 & " CONSTR_HOMOLOGY",constr_homology
99 refstr = index(controlcard,'REFSTR').gt.0
100 pdbref = index(controlcard,'PDBREF').gt.0
101 dyn_ss=(index(controlcard,'DYN_SS').gt.0)
102 C /06/28/2013 Adasko: dyn_ss is keyword allowing to break and create bond
103 C disulfide bond. Note that in conterary to dynamics this in
104 C CONTROLCARD. The bond are read in molread_zs.F
108 c------------------------------------------------------------------------------
109 subroutine read_efree(*)
111 C Read molecular data
115 include 'DIMENSIONS.ZSCOPT'
116 include 'DIMENSIONS.COMPAR'
117 include 'DIMENSIONS.FREE'
118 include 'COMMON.IOUNITS'
119 include 'COMMON.TIME1'
120 include 'COMMON.SBRIDGE'
121 include 'COMMON.CONTROL'
122 include 'COMMON.CHAIN'
123 include 'COMMON.HEADER'
125 include 'COMMON.FREE'
126 character*320 controlcard,ucase
127 integer iparm,ib,i,j,npars
139 call card_concat(controlcard,.true.)
140 call readi(controlcard,'NT',nT_h(iparm),1)
141 write (iout,*) "iparm",iparm," nt",nT_h(iparm)
143 if (nT_h(iparm).gt.MaxT_h) then
144 write (iout,*) "Error: parameter out of range: NT",nT_h(iparm),
148 replica(iparm)=index(controlcard,"REPLICA").gt.0
149 umbrella(iparm)=index(controlcard,"UMBRELLA").gt.0
150 if (umbrella(iparm) .and. homol_nset.gt.1) then
151 umbrella(iparm) = .false.
153 & "Replica in homology restraints weights UMBRELLA ignored,",iparm
155 read_iset(iparm)=index(controlcard,"READ_ISET").gt.0
156 write (iout,*) "nQ",nQ," nT",nT_h(iparm)," replica ",
157 & replica(iparm)," umbrella ",umbrella(iparm),
158 & " read_iset",read_iset(iparm)
161 call card_concat(controlcard,.true.)
162 call readi(controlcard,'NR',nR(ib,iparm),1)
163 if (umbrella(iparm)) then
166 nRR(ib,iparm)=nR(ib,iparm)
168 if (nR(ib,iparm).gt.MaxR) then
169 write (iout,*) "Error: parameter out of range: NR",
173 call reada(controlcard,'TEMP',beta_h(ib,iparm),298.0d0)
174 beta_h(ib,iparm)=1.0d0/(beta_h(ib,iparm)*1.987D-3)
175 call multreada(controlcard,'FI',f(1,ib,iparm),nR(ib,iparm),
178 call card_concat(controlcard,.true.)
179 call multreada(controlcard,'KH',KH(1,i,ib,iparm),nQ,
181 call multreada(controlcard,'Q0',Q0(1,i,ib,iparm),nQ,
186 write (iout,*) "ib",ib," beta_h",
187 & 1.0d0/(0.001987*beta_h(ib,iparm))
188 write (iout,*) "nR",nR(ib,iparm)
189 write (iout,*) "fi",(f(i,ib,iparm),i=1,nR(ib,iparm))
191 write (iout,*) "i",i," Kh",(Kh(j,i,ib,iparm),j=1,nQ),
192 & "q0",(q0(j,i,ib,iparm),j=1,nQ)
204 nR(ib,iparm)=nR(ib,1)
205 if (umbrella(iparm)) then
208 nRR(ib,iparm)=nR(ib,1)
210 beta_h(ib,iparm)=beta_h(ib,1)
212 f(i,ib,iparm)=f(i,ib,1)
214 KH(j,i,ib,iparm)=KH(j,i,ib,1)
215 Q0(j,i,ib,iparm)=Q0(j,i,ib,1)
218 replica(iparm)=replica(1)
219 umbrella(iparm)=umbrella(1)
220 read_iset(iparm)=read_iset(1)
228 c-----------------------------------------------------------------------------
229 subroutine read_protein_data(*)
232 include "DIMENSIONS.ZSCOPT"
233 include "DIMENSIONS.FREE"
236 integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
239 include "COMMON.CHAIN"
240 include "COMMON.IOUNITS"
241 include "COMMON.PROT"
242 include "COMMON.PROTFILES"
243 include "COMMON.NAMES"
244 include "COMMON.FREE"
245 include "COMMON.OBCINKA"
247 character*16000 controlcard
248 integer i,ii,ib,iR,iparm,ilen,iroof,nthr,npars
258 C Read names of files with conformation data.
259 if (replica(iparm)) then
265 do ii=1,nRR(ib,iparm)
266 write (iout,*) "Parameter set",iparm," temperature",ib,
269 call card_concat(controlcard,.true.)
270 write (iout,*) controlcard(:ilen(controlcard))
271 call readi(controlcard,"NFILE_BIN",nfile_bin(ii,ib,iparm),0)
272 call readi(controlcard,"NFILE_ASC",nfile_asc(ii,ib,iparm),0)
273 call readi(controlcard,"NFILE_CX",nfile_cx(ii,ib,iparm),0)
274 call readi(controlcard,"REC_START",rec_start(ii,ib,iparm),1)
275 call readi(controlcard,"REC_END",rec_end(ii,ib,iparm),
276 & maxstr*isampl(iparm)+rec_start(ii,ib,iparm)-1)
277 call reada(controlcard,"TIME_START",
278 & time_start_collect(ii,ib,iparm),0.0d0)
279 call reada(controlcard,"TIME_END",time_end_collect(ii,ib,iparm),
281 write (iout,*) "rec_start",rec_start(ii,ib,iparm),
282 & " rec_end",rec_end(ii,ib,iparm)
283 write (iout,*) "time_start",time_start_collect(ii,ib,iparm),
284 & " time_end",time_end_collect(ii,ib,iparm)
286 if (replica(iparm)) then
287 call readi(controlcard,"TOTRAJ",totraj(ii,iparm),1)
288 write (iout,*) "Number of trajectories",totraj(ii,iparm)
291 if (nfile_bin(ii,ib,iparm).lt.2
292 & .and. nfile_asc(ii,ib,iparm).eq.0
293 & .and. nfile_cx(ii,ib,iparm).eq.0) then
294 write (iout,*) "Error - no action specified!"
297 if (nfile_bin(ii,ib,iparm).gt.0) then
298 call card_concat(controlcard,.false.)
299 call split_string(controlcard,protfiles(1,1,ii,ib,iparm),
300 & maxfile_prot,nfile_bin(ii,ib,iparm))
302 write(iout,*)"nfile_bin",nfile_bin(ii,ib,iparm)
303 write(iout,*) (protfiles(i,1,ii,ib,iparm),
304 & i=1,nfile_bin(ii,ib,iparm))
307 if (nfile_asc(ii,ib,iparm).gt.0) then
308 call card_concat(controlcard,.false.)
309 call split_string(controlcard,protfiles(1,2,ii,ib,iparm),
310 & maxfile_prot,nfile_asc(ii,ib,iparm))
312 write(iout,*) "nfile_asc(ii,ib,iparm)",nfile_asc(ii,ib,iparm)
313 write(iout,*) (protfiles(i,2,ii,ib,iparm),
314 & i=1,nfile_asc(ii,ib,iparm))
316 else if (nfile_cx(ii,ib,iparm).gt.0) then
317 call card_concat(controlcard,.false.)
318 call split_string(controlcard,protfiles(1,2,ii,ib,iparm),
319 & maxfile_prot,nfile_cx(ii,ib,iparm))
321 write(iout,*) "nfile_cx(ii,ib,iparm)",nfile_cx(ii,ib,iparm)
322 write(iout,*) (protfiles(i,2,ii,ib,iparm),
323 & i=1,nfile_cx(ii,ib,iparm))
334 c-------------------------------------------------------------------------------
335 subroutine opentmp(islice,iunit,bprotfile_temp)
338 include "DIMENSIONS.ZSCOPT"
339 include "DIMENSIONS.FREE"
342 integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
345 include "COMMON.IOUNITS"
346 include "COMMON.PROTFILES"
347 include "COMMON.PROT"
348 include "COMMON.FREE"
349 character*64 bprotfile_temp
350 character*3 liczba,liczba2
357 write (liczba1,'(bz,i2.2)') islice
358 write (liczba,'(bz,i3.3)') me
360 c write (iout,*) "separate_parset ",separate_parset,
362 if (separate_parset) then
363 write (liczba2,'(bz,i3.3)') myparm
364 bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
365 & prefix(:ilen(prefix))//liczba//"_"//liczba2//".xbin.tmp"//liczba1
366 open (iunit,file=bprotfile_temp,status="unknown",
367 & form="unformatted",access="direct",recl=lenrec)
369 bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
370 & prefix(:ilen(prefix))//liczba//".xbin.tmp"//liczba1
371 open (iunit,file=bprotfile_temp,status="unknown",
372 & form="unformatted",access="direct",recl=lenrec)
375 bprotfile_temp = scratchdir(:ilen(scratchdir))//
376 & "/"//prefix(:ilen(prefix))//".xbin.tmp"//liczba1
377 open (iunit,file=bprotfile_temp,status="unknown",
378 & form="unformatted",access="direct",recl=lenrec)
380 c write (iout,*) "OpenTmp iunit",iunit," bprotfile_temp",
385 c-------------------------------------------------------------------------------
386 subroutine read_database(*)
389 include "DIMENSIONS.ZSCOPT"
390 include "DIMENSIONS.FREE"
393 integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
396 include "COMMON.CHAIN"
397 include "COMMON.IOUNITS"
398 include "COMMON.PROTFILES"
399 include "COMMON.NAMES"
402 include "COMMON.ENEPS"
403 include "COMMON.PROT"
404 include "COMMON.INTERACT"
405 include "COMMON.FREE"
406 include "COMMON.SBRIDGE"
407 include "COMMON.OBCINKA"
408 real*4 csingle(3,maxres2)
409 character*64 nazwa,bprotfile_temp
412 integer i,j,ii,jj(maxslice),k,kk(maxslice),l,
413 & ll(maxslice),mm(maxslice),if
414 integer nrec,nlines,iscor,iunit,islice
415 double precision energ
418 double precision rmsdev,energia(0:max_ene),efree,eini,temp
419 double precision prop(maxQ)
420 integer ntot_all(maxslice,0:maxprocs-1)
421 integer iparm,ib,iib,ir,nprop,nthr,npars
422 double precision etot,time
426 lenrec1=12*(nres+nct-nnt+1)+4*(2*nss+2)+24
427 lenrec2=12*(nres+nct-nnt+1)+4*(2*nss+2)+24+8*nQ
429 write (iout,*) "lenrec",lenrec," lenrec1",lenrec1,
439 write (iout,*) "nparmset",nparmset
447 if (replica(iparm)) then
454 do iR=1,nRR(ib,iparm)
456 write (iout,*) "iparm",iparm," ib",ib," iR",iR," nQ",nQ
462 IF (NFILE_BIN(iR,ib,iparm).GT.0) THEN
463 c Read conformations from binary DA files (one per batch) and write them to
464 c a binary DA scratchfile.
465 write (liczba,'(bz,i3.3)') me
466 do if=1,nfile_bin(iR,ib,iparm)
467 nazwa=protfiles(if,1,iR,ib,iparm)
468 & (:ilen(protfiles(if,1,iR,ib,iparm)))//".bx"
469 open (ientin,file=nazwa,status="old",form="unformatted",
470 & access="direct",recl=lenrec2,err=1111)
473 call opentmp(islice,ientout,bprotfile_temp)
474 call bxread(nazwa,ii,jj(islice),kk(islice),ll(islice),
475 & mm(islice),iR,ib,iparm)
482 IF (NFILE_ASC(iR,ib,iparm).GT.0) THEN
483 c Read conformations from multiple ASCII int files and write them to a binary
485 do if=1,nfile_asc(iR,ib,iparm)
486 nazwa=protfiles(if,2,iR,ib,iparm)
487 & (:ilen(protfiles(if,2,iR,ib,iparm)))//".x"
488 open(unit=ientin,file=nazwa,status='old',err=1111)
489 write(iout,*) "reading ",nazwa(:ilen(nazwa))
491 call xread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm)
494 IF (NFILE_CX(iR,ib,iparm).gt.0) THEN
495 c Read conformations from cx files and write them to a binary
497 do if=1,nfile_cx(iR,ib,iparm)
498 nazwa=protfiles(if,2,iR,ib,iparm)
499 & (:ilen(protfiles(if,2,iR,ib,iparm)))//".cx"
500 write(iout,*) "reading ",nazwa(:ilen(nazwa))
502 print *,"Calling cxread"
503 call cxread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm,
506 write (iout,*) "exit cxread"
512 stot(islice)=stot(islice)+jj(islice)
517 write (iout,*) "IPARM",iparm
520 if (nslice.eq.1) then
522 write (liczba,'(bz,i3.3)') me
523 bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
524 & prefix(:ilen(prefix))//liczba//".xbin.tmp"
526 bprotfile_temp = scratchdir(:ilen(scratchdir))//
527 & "/"//prefix(:ilen(prefix))//".xbin.tmp"
529 write(iout,*) mm(1)," conformations read",ll(1),
530 & " conformations written to ",
531 & bprotfile_temp(:ilen(bprotfile_temp))
534 write (liczba1,'(bz,i2.2)') islice
536 write (liczba,'(bz,i3.3)') me
537 bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
538 & prefix(:ilen(prefix))//liczba//".xbin.tmp"//liczba1
540 bprotfile_temp = scratchdir(:ilen(scratchdir))//
541 & "/"//prefix(:ilen(prefix))//".xbin.tmp"//liczba1
543 write(iout,*) mm(islice)," conformations read",ll(islice),
544 & " conformations written to ",
545 & bprotfile_temp(:ilen(bprotfile_temp))
550 c Check if everyone has the same number of conformations
551 call MPI_Allgather(stot(1),maxslice,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