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 read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
100 refstr = index(controlcard,'REFSTR').gt.0
101 pdbref = index(controlcard,'PDBREF').gt.0
102 dyn_ss=(index(controlcard,'DYN_SS').gt.0)
103 C /06/28/2013 Adasko: dyn_ss is keyword allowing to break and create bond
104 C disulfide bond. Note that in conterary to dynamics this in
105 C CONTROLCARD. The bond are read in molread_zs.F
109 c------------------------------------------------------------------------------
110 subroutine read_efree(*)
112 C Read molecular data
116 include 'DIMENSIONS.ZSCOPT'
117 include 'DIMENSIONS.COMPAR'
118 include 'DIMENSIONS.FREE'
119 include 'COMMON.IOUNITS'
120 include 'COMMON.TIME1'
121 include 'COMMON.SBRIDGE'
122 include 'COMMON.CONTROL'
123 include 'COMMON.CHAIN'
124 include 'COMMON.HEADER'
126 include 'COMMON.FREE'
127 character*320 controlcard,ucase
128 integer iparm,ib,i,j,npars
140 call card_concat(controlcard,.true.)
141 call readi(controlcard,'NT',nT_h(iparm),1)
142 write (iout,*) "iparm",iparm," nt",nT_h(iparm)
144 if (nT_h(iparm).gt.MaxT_h) then
145 write (iout,*) "Error: parameter out of range: NT",nT_h(iparm),
149 replica(iparm)=index(controlcard,"REPLICA").gt.0
150 umbrella(iparm)=index(controlcard,"UMBRELLA").gt.0
151 if (umbrella(iparm) .and. homol_nset.gt.1) then
152 umbrella(iparm) = .false.
154 & "Replica in homology restraints weights UMBRELLA ignored,",iparm
156 read_iset(iparm)=index(controlcard,"READ_ISET").gt.0
157 write (iout,*) "nQ",nQ," nT",nT_h(iparm)," replica ",
158 & replica(iparm)," umbrella ",umbrella(iparm),
159 & " read_iset",read_iset(iparm)
162 call card_concat(controlcard,.true.)
163 call readi(controlcard,'NR',nR(ib,iparm),1)
164 if (umbrella(iparm)) then
167 nRR(ib,iparm)=nR(ib,iparm)
169 if (nR(ib,iparm).gt.MaxR) then
170 write (iout,*) "Error: parameter out of range: NR",
174 call reada(controlcard,'TEMP',beta_h(ib,iparm),298.0d0)
175 beta_h(ib,iparm)=1.0d0/(beta_h(ib,iparm)*1.987D-3)
176 call multreada(controlcard,'FI',f(1,ib,iparm),nR(ib,iparm),
179 call card_concat(controlcard,.true.)
180 call multreada(controlcard,'KH',KH(1,i,ib,iparm),nQ,
182 call multreada(controlcard,'Q0',Q0(1,i,ib,iparm),nQ,
187 write (iout,*) "ib",ib," beta_h",
188 & 1.0d0/(0.001987*beta_h(ib,iparm))
189 write (iout,*) "nR",nR(ib,iparm)
190 write (iout,*) "fi",(f(i,ib,iparm),i=1,nR(ib,iparm))
192 write (iout,*) "i",i," Kh",(Kh(j,i,ib,iparm),j=1,nQ),
193 & "q0",(q0(j,i,ib,iparm),j=1,nQ)
198 write (iout,*) "HOMOL_NSET",homol_nset
199 if (homol_nset.gt.1) then
200 write (iout,*) "HOMOL_NSET: nT_h",nT_h(iparm)
202 nR(ib,iparm)=homol_nset
203 write (iout,*) "iparm",iparm," ib",ib," nR",nR(ib,iparm)
214 nR(ib,iparm)=nR(ib,1)
215 if (umbrella(iparm)) then
218 nRR(ib,iparm)=nR(ib,1)
220 beta_h(ib,iparm)=beta_h(ib,1)
222 f(i,ib,iparm)=f(i,ib,1)
224 KH(j,i,ib,iparm)=KH(j,i,ib,1)
225 Q0(j,i,ib,iparm)=Q0(j,i,ib,1)
228 replica(iparm)=replica(1)
229 umbrella(iparm)=umbrella(1)
230 read_iset(iparm)=read_iset(1)
238 c-----------------------------------------------------------------------------
239 subroutine read_protein_data(*)
242 include "DIMENSIONS.ZSCOPT"
243 include "DIMENSIONS.FREE"
246 integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
249 include "COMMON.CHAIN"
250 include "COMMON.IOUNITS"
251 include "COMMON.PROT"
252 include "COMMON.PROTFILES"
253 include "COMMON.NAMES"
254 include "COMMON.FREE"
255 include "COMMON.OBCINKA"
257 character*16000 controlcard
258 integer i,ii,ib,iR,iparm,ilen,iroof,nthr,npars
268 C Read names of files with conformation data.
269 if (replica(iparm)) then
275 do ii=1,nRR(ib,iparm)
276 write (iout,*) "Parameter set",iparm," temperature",ib,
279 call card_concat(controlcard,.true.)
280 write (iout,*) controlcard(:ilen(controlcard))
281 call readi(controlcard,"NFILE_BIN",nfile_bin(ii,ib,iparm),0)
282 call readi(controlcard,"NFILE_ASC",nfile_asc(ii,ib,iparm),0)
283 call readi(controlcard,"NFILE_CX",nfile_cx(ii,ib,iparm),0)
284 call readi(controlcard,"REC_START",rec_start(ii,ib,iparm),1)
285 call readi(controlcard,"REC_END",rec_end(ii,ib,iparm),
286 & maxstr*isampl(iparm)+rec_start(ii,ib,iparm)-1)
287 call reada(controlcard,"TIME_START",
288 & time_start_collect(ii,ib,iparm),0.0d0)
289 call reada(controlcard,"TIME_END",time_end_collect(ii,ib,iparm),
291 write (iout,*) "rec_start",rec_start(ii,ib,iparm),
292 & " rec_end",rec_end(ii,ib,iparm)
293 write (iout,*) "time_start",time_start_collect(ii,ib,iparm),
294 & " time_end",time_end_collect(ii,ib,iparm)
296 if (replica(iparm)) then
297 call readi(controlcard,"TOTRAJ",totraj(ii,iparm),1)
298 write (iout,*) "Number of trajectories",totraj(ii,iparm)
301 if (nfile_bin(ii,ib,iparm).lt.2
302 & .and. nfile_asc(ii,ib,iparm).eq.0
303 & .and. nfile_cx(ii,ib,iparm).eq.0) then
304 write (iout,*) "Error - no action specified!"
307 if (nfile_bin(ii,ib,iparm).gt.0) then
308 call card_concat(controlcard,.false.)
309 call split_string(controlcard,protfiles(1,1,ii,ib,iparm),
310 & maxfile_prot,nfile_bin(ii,ib,iparm))
312 write(iout,*)"nfile_bin",nfile_bin(ii,ib,iparm)
313 write(iout,*) (protfiles(i,1,ii,ib,iparm),
314 & i=1,nfile_bin(ii,ib,iparm))
317 if (nfile_asc(ii,ib,iparm).gt.0) then
318 call card_concat(controlcard,.false.)
319 call split_string(controlcard,protfiles(1,2,ii,ib,iparm),
320 & maxfile_prot,nfile_asc(ii,ib,iparm))
322 write(iout,*) "nfile_asc(ii,ib,iparm)",nfile_asc(ii,ib,iparm)
323 write(iout,*) (protfiles(i,2,ii,ib,iparm),
324 & i=1,nfile_asc(ii,ib,iparm))
326 else if (nfile_cx(ii,ib,iparm).gt.0) then
327 call card_concat(controlcard,.false.)
328 call split_string(controlcard,protfiles(1,2,ii,ib,iparm),
329 & maxfile_prot,nfile_cx(ii,ib,iparm))
331 write(iout,*) "nfile_cx(ii,ib,iparm)",nfile_cx(ii,ib,iparm)
332 write(iout,*) (protfiles(i,2,ii,ib,iparm),
333 & i=1,nfile_cx(ii,ib,iparm))
344 c-------------------------------------------------------------------------------
345 subroutine opentmp(islice,iunit,bprotfile_temp)
348 include "DIMENSIONS.ZSCOPT"
349 include "DIMENSIONS.FREE"
352 integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
355 include "COMMON.IOUNITS"
356 include "COMMON.PROTFILES"
357 include "COMMON.PROT"
358 include "COMMON.FREE"
359 character*128 bprotfile_temp
360 character*3 liczba,liczba2
367 write (liczba1,'(bz,i2.2)') islice
368 write (liczba,'(bz,i3.3)') me
370 c write (iout,*) "separate_parset ",separate_parset,
372 if (separate_parset) then
373 write (liczba2,'(bz,i3.3)') myparm
374 bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
375 & prefix(:ilen(prefix))//liczba//"_"//liczba2//".xbin.tmp"//liczba1
376 open (iunit,file=bprotfile_temp,status="unknown",
377 & form="unformatted",access="direct",recl=lenrec)
379 bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
380 & prefix(:ilen(prefix))//liczba//".xbin.tmp"//liczba1
381 open (iunit,file=bprotfile_temp,status="unknown",
382 & form="unformatted",access="direct",recl=lenrec)
385 bprotfile_temp = scratchdir(:ilen(scratchdir))//
386 & "/"//prefix(:ilen(prefix))//".xbin.tmp"//liczba1
387 open (iunit,file=bprotfile_temp,status="unknown",
388 & form="unformatted",access="direct",recl=lenrec)
390 c write (iout,*) "OpenTmp iunit",iunit," bprotfile_temp",
395 c-------------------------------------------------------------------------------
396 subroutine read_database(*)
399 include "DIMENSIONS.ZSCOPT"
400 include "DIMENSIONS.FREE"
403 integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
406 include "COMMON.CHAIN"
407 include "COMMON.IOUNITS"
408 include "COMMON.PROTFILES"
409 include "COMMON.NAMES"
412 include "COMMON.ENEPS"
413 include "COMMON.PROT"
414 include "COMMON.INTERACT"
415 include "COMMON.FREE"
416 include "COMMON.SBRIDGE"
417 include "COMMON.OBCINKA"
418 real*4 csingle(3,maxres2)
420 character*128 bprotfile_temp
423 integer i,j,ii,jj(maxslice),k,kk(maxslice),l,
424 & ll(maxslice),mm(maxslice),if
425 integer nrec,nlines,iscor,iunit,islice
426 double precision energ
429 double precision rmsdev,energia(0:max_ene),efree,eini,temp
430 double precision prop(maxQ)
431 integer ntot_all(maxslice,0:maxprocs-1)
432 integer iparm,ib,iib,ir,nprop,nthr,npars
433 double precision etot,time
437 lenrec1=12*(nres+nct-nnt+1)+4*(2*nss+2)+24
438 lenrec2=12*(nres+nct-nnt+1)+4*(2*nss+2)+24+8*nQ
440 write (iout,*) "lenrec",lenrec," lenrec1",lenrec1,
450 write (iout,*) "nparmset",nparmset
458 if (replica(iparm)) then
465 do iR=1,nRR(ib,iparm)
467 write (iout,*) "iparm",iparm," ib",ib," iR",iR," nQ",nQ
473 IF (NFILE_BIN(iR,ib,iparm).GT.0) THEN
474 c Read conformations from binary DA files (one per batch) and write them to
475 c a binary DA scratchfile.
476 write (liczba,'(bz,i3.3)') me
477 do if=1,nfile_bin(iR,ib,iparm)
478 nazwa=protfiles(if,1,iR,ib,iparm)
479 & (:ilen(protfiles(if,1,iR,ib,iparm)))//".bx"
480 open (ientin,file=nazwa,status="old",form="unformatted",
481 & access="direct",recl=lenrec2,err=1111)
484 call opentmp(islice,ientout,bprotfile_temp)
485 call bxread(nazwa,ii,jj(islice),kk(islice),ll(islice),
486 & mm(islice),iR,ib,iparm)
493 IF (NFILE_ASC(iR,ib,iparm).GT.0) THEN
494 c Read conformations from multiple ASCII int files and write them to a binary
496 do if=1,nfile_asc(iR,ib,iparm)
497 nazwa=protfiles(if,2,iR,ib,iparm)
498 & (:ilen(protfiles(if,2,iR,ib,iparm)))//".x"
499 open(unit=ientin,file=nazwa,status='old',err=1111)
500 write(iout,*) "reading ",nazwa(:ilen(nazwa))
502 call xread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm)
505 IF (NFILE_CX(iR,ib,iparm).gt.0) THEN
506 c Read conformations from cx files and write them to a binary
508 do if=1,nfile_cx(iR,ib,iparm)
509 nazwa=protfiles(if,2,iR,ib,iparm)
510 & (:ilen(protfiles(if,2,iR,ib,iparm)))//".cx"
511 write(iout,*) "reading ",nazwa(:ilen(nazwa))
513 print *,"Calling cxread"
514 call cxread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm,
517 write (iout,*) "exit cxread"
523 stot(islice)=stot(islice)+jj(islice)
528 write (iout,*) "IPARM",iparm
531 if (nslice.eq.1) then
533 write (liczba,'(bz,i3.3)') me
534 bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
535 & prefix(:ilen(prefix))//liczba//".xbin.tmp"
537 bprotfile_temp = scratchdir(:ilen(scratchdir))//
538 & "/"//prefix(:ilen(prefix))//".xbin.tmp"
540 write(iout,*) mm(1)," conformations read",ll(1),
541 & " conformations written to ",
542 & bprotfile_temp(:ilen(bprotfile_temp))
545 write (liczba1,'(bz,i2.2)') islice
547 write (liczba,'(bz,i3.3)') me
548 bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
549 & prefix(:ilen(prefix))//liczba//".xbin.tmp"//liczba1
551 bprotfile_temp = scratchdir(:ilen(scratchdir))//
552 & "/"//prefix(:ilen(prefix))//".xbin.tmp"//liczba1
554 write(iout,*) mm(islice)," conformations read",ll(islice),
555 & " conformations written to ",
556 & bprotfile_temp(:ilen(bprotfile_temp))
561 c Check if everyone has the same number of conformations
562 call MPI_Allgather(stot(1),maxslice,MPI_INTEGER,
563 & ntot_all(1,0),maxslice,MPI_INTEGER,MPI_Comm_World,IERROR)
568 if (stot(islice).ne.ntot_all(islice,i)) then
569 write (iout,*) "Number of conformations at processor",i,
570 & " differs from that at processor",me,
571 & stot(islice),ntot_all(islice,i)," slice",islice
579 write (iout,*) "Numbers of conformations read by processors"
582 write (iout,'(8i10)') i,(ntot_all(islice,i),islice=1,nslice)
584 write (iout,*) "Calculation terminated."
589 ntot(islice)=stot(islice)
593 1111 write(iout,*) "Error opening coordinate file ",nazwa(:ilen(nazwa))
597 c------------------------------------------------------------------------------
598 subroutine card_concat(card,to_upper)
600 include 'DIMENSIONS.ZSCOPT'
601 include "COMMON.IOUNITS"
603 character*80 karta,ucase
607 read (inp,'(a)') karta
608 if (to_upper) karta=ucase(karta)
610 do while (karta(80:80).eq.'&')
611 card=card(:ilen(card)+1)//karta(:79)
612 read (inp,'(a)') karta
613 if (to_upper) karta=ucase(karta)
615 card=card(:ilen(card)+1)//karta
618 c------------------------------------------------------------------------------
619 subroutine readi(rekord,lancuch,wartosc,default)
621 character*(*) rekord,lancuch
622 integer wartosc,default
625 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
630 iread=iread+ilen(lancuch)+1
631 read (rekord(iread:),*) wartosc
634 c----------------------------------------------------------------------------
635 subroutine reada(rekord,lancuch,wartosc,default)
637 character*(*) rekord,lancuch
639 double precision wartosc,default
642 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
647 iread=iread+ilen(lancuch)+1
648 read (rekord(iread:),*) wartosc
651 c----------------------------------------------------------------------------
652 subroutine multreadi(rekord,lancuch,tablica,dim,default)
655 integer tablica(dim),default
656 character*(*) rekord,lancuch
663 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
664 if (iread.eq.0) return
665 iread=iread+ilen(lancuch)+1
666 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
669 c----------------------------------------------------------------------------
670 subroutine multreada(rekord,lancuch,tablica,dim,default)
673 double precision tablica(dim),default
674 character*(*) rekord,lancuch
681 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
682 if (iread.eq.0) return
683 iread=iread+ilen(lancuch)+1
684 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
687 c----------------------------------------------------------------------------
688 subroutine reads(rekord,lancuch,wartosc,default)
690 character*(*) rekord,lancuch,wartosc,default
692 integer ilen,lenlan,lenrec,iread,ireade
698 iread=index(rekord,lancuch(:lenlan)//"=")
699 c print *,"rekord",rekord," lancuch",lancuch
700 c print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
706 c print *,"iread",iread
707 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
708 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
710 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
712 c print *,"iread",iread
713 if (iread.gt.lenrec) then
718 c print *,"ireade",ireade
719 do while (ireade.lt.lenrec .and.
720 & .not.iblnk(rekord(ireade:ireade)))
723 wartosc=rekord(iread:ireade)
726 c----------------------------------------------------------------------------
727 subroutine multreads(rekord,lancuch,tablica,dim,default)
730 character*(*) rekord,lancuch,tablica(dim),default
732 integer ilen,lenlan,lenrec,iread,ireade
741 iread=index(rekord,lancuch(:lenlan)//"=")
742 c print *,"rekord",rekord," lancuch",lancuch
743 c print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
744 if (iread.eq.0) return
747 c print *,"iread",iread
748 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
749 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
751 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
753 c print *,"iread",iread
754 if (iread.gt.lenrec) return
756 c print *,"ireade",ireade
757 do while (ireade.lt.lenrec .and.
758 & .not.iblnk(rekord(ireade:ireade)))
761 tablica(i)=rekord(iread:ireade)
765 c----------------------------------------------------------------------------
766 subroutine split_string(rekord,tablica,dim,nsub)
768 integer dim,nsub,i,ii,ll,kk
769 character*(*) tablica(dim)
780 C Find the start of term name
782 do while (ii.le.ll .and. rekord(ii:ii).eq." ")
785 C Parse the name into TABLICA(i) until blank found
786 do while (ii.le.ll .and. rekord(ii:ii).ne." ")
788 tablica(i)(kk:kk)=rekord(ii:ii)
791 if (kk.gt.0) nsub=nsub+1
796 c--------------------------------------------------------------------------------
797 integer function iroof(n,m)
799 if (ii*m .lt. n) ii=ii+1