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 refstr = index(controlcard,'REFSTR').gt.0
97 pdbref = index(controlcard,'PDBREF').gt.0
98 dyn_ss=(index(controlcard,'DYN_SS').gt.0)
99 C /06/28/2013 Adasko: dyn_ss is keyword allowing to break and create bond
100 C disulfide bond. Note that in conterary to dynamics this in
101 C CONTROLCARD. The bond are read in molread_zs.F
105 c------------------------------------------------------------------------------
106 subroutine read_efree(*)
108 C Read molecular data
112 include 'DIMENSIONS.ZSCOPT'
113 include 'DIMENSIONS.COMPAR'
114 include 'DIMENSIONS.FREE'
115 include 'COMMON.IOUNITS'
116 include 'COMMON.TIME1'
117 include 'COMMON.SBRIDGE'
118 include 'COMMON.CONTROL'
119 include 'COMMON.CHAIN'
120 include 'COMMON.HEADER'
122 include 'COMMON.FREE'
123 character*320 controlcard,ucase
124 integer iparm,ib,i,j,npars
136 call card_concat(controlcard,.true.)
137 call readi(controlcard,'NT',nT_h(iparm),1)
138 write (iout,*) "iparm",iparm," nt",nT_h(iparm)
140 if (nT_h(iparm).gt.MaxT_h) then
141 write (iout,*) "Error: parameter out of range: NT",nT_h(iparm),
145 replica(iparm)=index(controlcard,"REPLICA").gt.0
146 umbrella(iparm)=index(controlcard,"UMBRELLA").gt.0
147 read_iset(iparm)=index(controlcard,"READ_ISET").gt.0
148 write (iout,*) "nQ",nQ," nT",nT_h(iparm)," replica ",
149 & replica(iparm)," umbrella ",umbrella(iparm),
150 & " read_iset",read_iset(iparm)
153 call card_concat(controlcard,.true.)
154 call readi(controlcard,'NR',nR(ib,iparm),1)
155 if (umbrella(iparm)) then
158 nRR(ib,iparm)=nR(ib,iparm)
160 if (nR(ib,iparm).gt.MaxR) then
161 write (iout,*) "Error: parameter out of range: NR",
165 call reada(controlcard,'TEMP',beta_h(ib,iparm),298.0d0)
166 beta_h(ib,iparm)=1.0d0/(beta_h(ib,iparm)*1.987D-3)
167 call multreada(controlcard,'FI',f(1,ib,iparm),nR(ib,iparm),
170 call card_concat(controlcard,.true.)
171 call multreada(controlcard,'KH',KH(1,i,ib,iparm),nQ,
173 call multreada(controlcard,'Q0',Q0(1,i,ib,iparm),nQ,
178 write (iout,*) "ib",ib," beta_h",
179 & 1.0d0/(0.001987*beta_h(ib,iparm))
180 write (iout,*) "nR",nR(ib,iparm)
181 write (iout,*) "fi",(f(i,ib,iparm),i=1,nR(ib,iparm))
183 write (iout,*) "i",i," Kh",(Kh(j,i,ib,iparm),j=1,nQ),
184 & "q0",(q0(j,i,ib,iparm),j=1,nQ)
196 nR(ib,iparm)=nR(ib,1)
197 if (umbrella(iparm)) then
200 nRR(ib,iparm)=nR(ib,1)
202 beta_h(ib,iparm)=beta_h(ib,1)
204 f(i,ib,iparm)=f(i,ib,1)
206 KH(j,i,ib,iparm)=KH(j,i,ib,1)
207 Q0(j,i,ib,iparm)=Q0(j,i,ib,1)
210 replica(iparm)=replica(1)
211 umbrella(iparm)=umbrella(1)
212 read_iset(iparm)=read_iset(1)
220 c-----------------------------------------------------------------------------
221 subroutine read_protein_data(*)
224 include "DIMENSIONS.ZSCOPT"
225 include "DIMENSIONS.FREE"
228 integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
231 include "COMMON.CHAIN"
232 include "COMMON.IOUNITS"
233 include "COMMON.PROT"
234 include "COMMON.PROTFILES"
235 include "COMMON.NAMES"
236 include "COMMON.FREE"
237 include "COMMON.OBCINKA"
239 character*16000 controlcard
240 integer i,ii,ib,iR,iparm,ilen,iroof,nthr,npars
250 C Read names of files with conformation data.
251 if (replica(iparm)) then
257 do ii=1,nRR(ib,iparm)
258 write (iout,*) "Parameter set",iparm," temperature",ib,
261 call card_concat(controlcard,.true.)
262 write (iout,*) controlcard(:ilen(controlcard))
263 call readi(controlcard,"NFILE_BIN",nfile_bin(ii,ib,iparm),0)
264 call readi(controlcard,"NFILE_ASC",nfile_asc(ii,ib,iparm),0)
265 call readi(controlcard,"NFILE_CX",nfile_cx(ii,ib,iparm),0)
266 call readi(controlcard,"REC_START",rec_start(ii,ib,iparm),1)
267 call readi(controlcard,"REC_END",rec_end(ii,ib,iparm),
268 & maxstr*isampl(iparm)+rec_start(ii,ib,iparm)-1)
269 call reada(controlcard,"TIME_START",
270 & time_start_collect(ii,ib,iparm),0.0d0)
271 call reada(controlcard,"TIME_END",time_end_collect(ii,ib,iparm),
273 write (iout,*) "rec_start",rec_start(ii,ib,iparm),
274 & " rec_end",rec_end(ii,ib,iparm)
275 write (iout,*) "time_start",time_start_collect(ii,ib,iparm),
276 & " time_end",time_end_collect(ii,ib,iparm)
278 if (replica(iparm)) then
279 call readi(controlcard,"TOTRAJ",totraj(ii,iparm),1)
280 write (iout,*) "Number of trajectories",totraj(ii,iparm)
283 if (nfile_bin(ii,ib,iparm).lt.2
284 & .and. nfile_asc(ii,ib,iparm).eq.0
285 & .and. nfile_cx(ii,ib,iparm).eq.0) then
286 write (iout,*) "Error - no action specified!"
289 if (nfile_bin(ii,ib,iparm).gt.0) then
290 call card_concat(controlcard,.false.)
291 call split_string(controlcard,protfiles(1,1,ii,ib,iparm),
292 & maxfile_prot,nfile_bin(ii,ib,iparm))
294 write(iout,*)"nfile_bin",nfile_bin(ii,ib,iparm)
295 write(iout,*) (protfiles(i,1,ii,ib,iparm),
296 & i=1,nfile_bin(ii,ib,iparm))
299 if (nfile_asc(ii,ib,iparm).gt.0) then
300 call card_concat(controlcard,.false.)
301 call split_string(controlcard,protfiles(1,2,ii,ib,iparm),
302 & maxfile_prot,nfile_asc(ii,ib,iparm))
304 write(iout,*) "nfile_asc(ii,ib,iparm)",nfile_asc(ii,ib,iparm)
305 write(iout,*) (protfiles(i,2,ii,ib,iparm),
306 & i=1,nfile_asc(ii,ib,iparm))
308 else if (nfile_cx(ii,ib,iparm).gt.0) then
309 call card_concat(controlcard,.false.)
310 call split_string(controlcard,protfiles(1,2,ii,ib,iparm),
311 & maxfile_prot,nfile_cx(ii,ib,iparm))
313 write(iout,*) "nfile_cx(ii,ib,iparm)",nfile_cx(ii,ib,iparm)
314 write(iout,*) (protfiles(i,2,ii,ib,iparm),
315 & i=1,nfile_cx(ii,ib,iparm))
326 c-------------------------------------------------------------------------------
327 subroutine opentmp(islice,iunit,bprotfile_temp)
330 include "DIMENSIONS.ZSCOPT"
331 include "DIMENSIONS.FREE"
334 integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
337 include "COMMON.IOUNITS"
338 include "COMMON.PROTFILES"
339 include "COMMON.PROT"
340 include "COMMON.FREE"
341 character*64 bprotfile_temp
342 character*3 liczba,liczba2
349 write (liczba1,'(bz,i2.2)') islice
350 write (liczba,'(bz,i3.3)') me
352 c write (iout,*) "separate_parset ",separate_parset,
354 if (separate_parset) then
355 write (liczba2,'(bz,i3.3)') myparm
356 bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
357 & prefix(:ilen(prefix))//liczba//"_"//liczba2//".xbin.tmp"//liczba1
358 open (iunit,file=bprotfile_temp,status="unknown",
359 & form="unformatted",access="direct",recl=lenrec)
361 bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
362 & prefix(:ilen(prefix))//liczba//".xbin.tmp"//liczba1
363 open (iunit,file=bprotfile_temp,status="unknown",
364 & form="unformatted",access="direct",recl=lenrec)
367 bprotfile_temp = scratchdir(:ilen(scratchdir))//
368 & "/"//prefix(:ilen(prefix))//".xbin.tmp"//liczba1
369 open (iunit,file=bprotfile_temp,status="unknown",
370 & form="unformatted",access="direct",recl=lenrec)
372 c write (iout,*) "OpenTmp iunit",iunit," bprotfile_temp",
377 c-------------------------------------------------------------------------------
378 subroutine read_database(*)
381 include "DIMENSIONS.ZSCOPT"
382 include "DIMENSIONS.FREE"
385 integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
388 include "COMMON.CHAIN"
389 include "COMMON.IOUNITS"
390 include "COMMON.PROTFILES"
391 include "COMMON.NAMES"
394 include "COMMON.ENEPS"
395 include "COMMON.PROT"
396 include "COMMON.INTERACT"
397 include "COMMON.FREE"
398 include "COMMON.SBRIDGE"
399 include "COMMON.OBCINKA"
400 real*4 csingle(3,maxres2)
401 character*64 nazwa,bprotfile_temp
404 integer i,j,ii,jj(maxslice),k,kk(maxslice),l,
405 & ll(maxslice),mm(maxslice),if
406 integer nrec,nlines,iscor,iunit,islice
407 double precision energ
410 double precision rmsdev,energia(0:max_ene),efree,eini,temp
411 double precision prop(maxQ)
412 integer ntot_all(maxslice,0:maxprocs-1)
413 integer iparm,ib,iib,ir,nprop,nthr,npars
414 double precision etot,time
418 lenrec1=12*(nres+nct-nnt+1)+4*(2*nss+2)+24
419 lenrec2=12*(nres+nct-nnt+1)+4*(2*nss+2)+24+8*nQ
421 write (iout,*) "lenrec",lenrec," lenrec1",lenrec1,
431 write (iout,*) "nparmset",nparmset
439 if (replica(iparm)) then
446 do iR=1,nRR(ib,iparm)
448 write (iout,*) "iparm",iparm," ib",ib," iR",iR," nQ",nQ
454 IF (NFILE_BIN(iR,ib,iparm).GT.0) THEN
455 c Read conformations from binary DA files (one per batch) and write them to
456 c a binary DA scratchfile.
457 write (liczba,'(bz,i3.3)') me
458 do if=1,nfile_bin(iR,ib,iparm)
459 nazwa=protfiles(if,1,iR,ib,iparm)
460 & (:ilen(protfiles(if,1,iR,ib,iparm)))//".bx"
461 open (ientin,file=nazwa,status="old",form="unformatted",
462 & access="direct",recl=lenrec2,err=1111)
465 call opentmp(islice,ientout,bprotfile_temp)
466 call bxread(nazwa,ii,jj(islice),kk(islice),ll(islice),
467 & mm(islice),iR,ib,iparm)
474 IF (NFILE_ASC(iR,ib,iparm).GT.0) THEN
475 c Read conformations from multiple ASCII int files and write them to a binary
477 do if=1,nfile_asc(iR,ib,iparm)
478 nazwa=protfiles(if,2,iR,ib,iparm)
479 & (:ilen(protfiles(if,2,iR,ib,iparm)))//".x"
480 open(unit=ientin,file=nazwa,status='old',err=1111)
481 write(iout,*) "reading ",nazwa(:ilen(nazwa))
483 call xread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm)
486 IF (NFILE_CX(iR,ib,iparm).gt.0) THEN
487 c Read conformations from cx files and write them to a binary
489 do if=1,nfile_cx(iR,ib,iparm)
490 nazwa=protfiles(if,2,iR,ib,iparm)
491 & (:ilen(protfiles(if,2,iR,ib,iparm)))//".cx"
492 write(iout,*) "reading ",nazwa(:ilen(nazwa))
494 print *,"Calling cxread"
495 call cxread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm,
498 write (iout,*) "exit cxread"
504 stot(islice)=stot(islice)+jj(islice)
509 write (iout,*) "IPARM",iparm
512 if (nslice.eq.1) then
514 write (liczba,'(bz,i3.3)') me
515 bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
516 & prefix(:ilen(prefix))//liczba//".xbin.tmp"
518 bprotfile_temp = scratchdir(:ilen(scratchdir))//
519 & "/"//prefix(:ilen(prefix))//".xbin.tmp"
521 write(iout,*) mm(1)," conformations read",ll(1),
522 & " conformations written to ",
523 & bprotfile_temp(:ilen(bprotfile_temp))
526 write (liczba1,'(bz,i2.2)') islice
528 write (liczba,'(bz,i3.3)') me
529 bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
530 & prefix(:ilen(prefix))//liczba//".xbin.tmp"//liczba1
532 bprotfile_temp = scratchdir(:ilen(scratchdir))//
533 & "/"//prefix(:ilen(prefix))//".xbin.tmp"//liczba1
535 write(iout,*) mm(islice)," conformations read",ll(islice),
536 & " conformations written to ",
537 & bprotfile_temp(:ilen(bprotfile_temp))
542 c Check if everyone has the same number of conformations
543 call MPI_Allgather(stot(1),maxslice,MPI_INTEGER,
544 & ntot_all(1,0),maxslice,MPI_INTEGER,MPI_Comm_World,IERROR)
549 if (stot(islice).ne.ntot_all(islice,i)) then
550 write (iout,*) "Number of conformations at processor",i,
551 & " differs from that at processor",me,
552 & stot(islice),ntot_all(islice,i)," slice",islice
560 write (iout,*) "Numbers of conformations read by processors"
563 write (iout,'(8i10)') i,(ntot_all(islice,i),islice=1,nslice)
565 write (iout,*) "Calculation terminated."
570 ntot(islice)=stot(islice)
574 1111 write(iout,*) "Error opening coordinate file ",nazwa(:ilen(nazwa))
578 c------------------------------------------------------------------------------
579 subroutine card_concat(card,to_upper)
581 include 'DIMENSIONS.ZSCOPT'
582 include "COMMON.IOUNITS"
584 character*80 karta,ucase
588 read (inp,'(a)') karta
589 if (to_upper) karta=ucase(karta)
591 do while (karta(80:80).eq.'&')
592 card=card(:ilen(card)+1)//karta(:79)
593 read (inp,'(a)') karta
594 if (to_upper) karta=ucase(karta)
596 card=card(:ilen(card)+1)//karta
599 c------------------------------------------------------------------------------
600 subroutine readi(rekord,lancuch,wartosc,default)
602 character*(*) rekord,lancuch
603 integer wartosc,default
606 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
611 iread=iread+ilen(lancuch)+1
612 read (rekord(iread:),*) wartosc
615 c----------------------------------------------------------------------------
616 subroutine reada(rekord,lancuch,wartosc,default)
618 character*(*) rekord,lancuch
620 double precision wartosc,default
623 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
628 iread=iread+ilen(lancuch)+1
629 read (rekord(iread:),*) wartosc
632 c----------------------------------------------------------------------------
633 subroutine multreadi(rekord,lancuch,tablica,dim,default)
636 integer tablica(dim),default
637 character*(*) rekord,lancuch
644 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
645 if (iread.eq.0) return
646 iread=iread+ilen(lancuch)+1
647 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
650 c----------------------------------------------------------------------------
651 subroutine multreada(rekord,lancuch,tablica,dim,default)
654 double precision tablica(dim),default
655 character*(*) rekord,lancuch
662 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
663 if (iread.eq.0) return
664 iread=iread+ilen(lancuch)+1
665 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
668 c----------------------------------------------------------------------------
669 subroutine reads(rekord,lancuch,wartosc,default)
671 character*(*) rekord,lancuch,wartosc,default
673 integer ilen,lenlan,lenrec,iread,ireade
679 iread=index(rekord,lancuch(:lenlan)//"=")
680 c print *,"rekord",rekord," lancuch",lancuch
681 c print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
687 c print *,"iread",iread
688 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
689 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
691 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
693 c print *,"iread",iread
694 if (iread.gt.lenrec) then
699 c print *,"ireade",ireade
700 do while (ireade.lt.lenrec .and.
701 & .not.iblnk(rekord(ireade:ireade)))
704 wartosc=rekord(iread:ireade)
707 c----------------------------------------------------------------------------
708 subroutine multreads(rekord,lancuch,tablica,dim,default)
711 character*(*) rekord,lancuch,tablica(dim),default
713 integer ilen,lenlan,lenrec,iread,ireade
722 iread=index(rekord,lancuch(:lenlan)//"=")
723 c print *,"rekord",rekord," lancuch",lancuch
724 c print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
725 if (iread.eq.0) return
728 c print *,"iread",iread
729 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
730 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
732 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
734 c print *,"iread",iread
735 if (iread.gt.lenrec) return
737 c print *,"ireade",ireade
738 do while (ireade.lt.lenrec .and.
739 & .not.iblnk(rekord(ireade:ireade)))
742 tablica(i)=rekord(iread:ireade)
746 c----------------------------------------------------------------------------
747 subroutine split_string(rekord,tablica,dim,nsub)
749 integer dim,nsub,i,ii,ll,kk
750 character*(*) tablica(dim)
761 C Find the start of term name
763 do while (ii.le.ll .and. rekord(ii:ii).eq." ")
766 C Parse the name into TABLICA(i) until blank found
767 do while (ii.le.ll .and. rekord(ii:ii).ne." ")
769 tablica(i)(kk:kk)=rekord(ii:ii)
772 if (kk.gt.0) nsub=nsub+1
777 c--------------------------------------------------------------------------------
778 integer function iroof(n,m)
780 if (ii*m .lt. n) ii=ii+1