Merge branch 'lipid' into AFM
[unres.git] / source / wham / src-M / readrtns.F
1       subroutine read_general_data(*)
2       implicit none
3       include "DIMENSIONS"
4       include "DIMENSIONS.ZSCOPT"
5       include "DIMENSIONS.FREE"
6       include "COMMON.TORSION"
7       include "COMMON.INTERACT"
8       include "COMMON.IOUNITS"
9       include "COMMON.TIME1"
10       include "COMMON.PROT"
11       include "COMMON.PROTFILES"
12       include "COMMON.CHAIN"
13       include "COMMON.NAMES"
14       include "COMMON.FFIELD"
15       include "COMMON.ENEPS"
16       include "COMMON.WEIGHTS"
17       include "COMMON.FREE"
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
25       integer ilen
26       external ilen
27       character*16 ucase
28       character*16 key
29       external ucase
30
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,
35      &    max_ene
36         return1
37       endif
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",
44      &    nparmset, Max_Parm
45         return1
46       endif
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)
55       call flush(iout)
56       if (nslice.gt.MaxSlice) then
57         write (iout,*) "Error: parameter out of range: NSLICE",nslice,
58      &    MaxSlice
59         return1
60       endif
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)
65       if (nQ.gt.MaxQ) then
66         write (iout,*) "Error: parameter out of range: NQ",nq,
67      &    maxq
68         return1
69       endif
70       indpdb=0
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 reada(controlcard,"LIPTHICK",lipthick,0.0d0)
86       call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
87       if (lipthick.gt.0.0d0) then
88        bordliptop=(boxzsize+lipthick)/2.0
89        bordlipbot=bordliptop-lipthick
90 C      endif
91       if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
92      & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
93       buflipbot=bordlipbot+lipbufthick
94       bufliptop=bordliptop-lipbufthick
95       if ((lipbufthick*2.0d0).gt.lipthick)
96      &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
97       endif
98       write(iout,*) "bordliptop=",bordliptop
99       write(iout,*) "bordlipbot=",bordlipbot
100       write(iout,*) "bufliptop=",bufliptop
101       write(iout,*) "buflipbot=",buflipbot
102       call readi(controlcard,'SYM',symetr,1)
103       write (iout,*) "DISTCHAINMAX",distchainmax
104       write (iout,*) "delta",delta
105       write (iout,*) "einicheck",einicheck
106       write (iout,*) "rescale_mode",rescale_mode
107       call flush(iout)
108       bxfile=index(controlcard,"BXFILE").gt.0
109       cxfile=index(controlcard,"CXFILE").gt.0
110       if (nslice .eq. 1 .and. .not.bxfile .and. .not.cxfile)
111      & bxfile=.true.
112       histfile=index(controlcard,"HISTFILE").gt.0
113       histout=index(controlcard,"HISTOUT").gt.0
114       entfile=index(controlcard,"ENTFILE").gt.0
115       zscfile=index(controlcard,"ZSCFILE").gt.0
116       with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
117       write (iout,*) "with_dihed_constr ",with_dihed_constr
118       with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
119       write (iout,*) "with_theta_constr ",with_theta_constr
120       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
121       dyn_ss=index(controlcard,"DYN_SS").gt.0
122       return
123       end
124 c------------------------------------------------------------------------------
125       subroutine read_efree(*)
126 C
127 C Read molecular data
128 C
129       implicit none
130       include 'DIMENSIONS'
131       include 'DIMENSIONS.ZSCOPT'
132       include 'DIMENSIONS.COMPAR'
133       include 'DIMENSIONS.FREE'
134       include 'COMMON.IOUNITS'
135       include 'COMMON.TIME1'
136       include 'COMMON.SBRIDGE'
137       include 'COMMON.CONTROL'
138       include 'COMMON.CHAIN'
139       include 'COMMON.HEADER'
140       include 'COMMON.GEO'
141       include 'COMMON.FREE'
142       character*320 controlcard,ucase
143       integer iparm,ib,i,j,npars
144       integer ilen
145       external ilen
146      
147       if (hamil_rep) then
148         npars=1
149       else
150         npars=nParmSet
151       endif
152
153       do iparm=1,npars
154
155       call card_concat(controlcard,.true.)
156       call readi(controlcard,'NT',nT_h(iparm),1)
157       write (iout,*) "iparm",iparm," nt",nT_h(iparm)
158       call flush(iout)
159       if (nT_h(iparm).gt.MaxT_h) then
160         write (iout,*)  "Error: parameter out of range: NT",nT_h(iparm),
161      &    MaxT_h
162         return1
163       endif
164       replica(iparm)=index(controlcard,"REPLICA").gt.0
165       umbrella(iparm)=index(controlcard,"UMBRELLA").gt.0
166       read_iset(iparm)=index(controlcard,"READ_ISET").gt.0
167       write (iout,*) "nQ",nQ," nT",nT_h(iparm)," replica ",
168      &  replica(iparm)," umbrella ",umbrella(iparm),
169      &  " read_iset",read_iset(iparm)
170       call flush(iout)
171       do ib=1,nT_h(iparm)
172         call card_concat(controlcard,.true.)
173         call readi(controlcard,'NR',nR(ib,iparm),1)
174         if (umbrella(iparm)) then
175           nRR(ib,iparm)=1
176         else
177           nRR(ib,iparm)=nR(ib,iparm)
178         endif
179         if (nR(ib,iparm).gt.MaxR) then
180           write (iout,*)  "Error: parameter out of range: NR",
181      &      nR(ib,iparm),MaxR
182         return1
183         endif
184         call reada(controlcard,'TEMP',beta_h(ib,iparm),298.0d0)
185         beta_h(ib,iparm)=1.0d0/(beta_h(ib,iparm)*1.987D-3)
186         call multreada(controlcard,'FI',f(1,ib,iparm),nR(ib,iparm),
187      &    0.0d0)
188         do i=1,nR(ib,iparm)
189           call card_concat(controlcard,.true.)
190           call multreada(controlcard,'KH',KH(1,i,ib,iparm),nQ,
191      &      100.0d0)
192           call multreada(controlcard,'Q0',Q0(1,i,ib,iparm),nQ,
193      &      0.0d0)
194         enddo
195       enddo
196       do ib=1,nT_h(iparm)
197         write (iout,*) "ib",ib," beta_h",
198      &    1.0d0/(0.001987*beta_h(ib,iparm))
199         write (iout,*) "nR",nR(ib,iparm)
200         write (iout,*) "fi",(f(i,ib,iparm),i=1,nR(ib,iparm))
201         do i=1,nR(ib,iparm)
202           write (iout,*) "i",i," Kh",(Kh(j,i,ib,iparm),j=1,nQ),
203      &      "q0",(q0(j,i,ib,iparm),j=1,nQ)
204         enddo
205         call flush(iout)
206       enddo
207
208       enddo
209
210       if (hamil_rep) then
211
212        do iparm=2,nParmSet
213           nT_h(iparm)=nT_h(1)
214          do ib=1,nT_h(iparm)
215            nR(ib,iparm)=nR(ib,1)
216            if (umbrella(iparm)) then
217              nRR(ib,iparm)=1
218            else
219              nRR(ib,iparm)=nR(ib,1)
220            endif
221            beta_h(ib,iparm)=beta_h(ib,1)
222            do i=1,nR(ib,iparm)
223              f(i,ib,iparm)=f(i,ib,1)
224              do j=1,nQ
225                KH(j,i,ib,iparm)=KH(j,i,ib,1) 
226                Q0(j,i,ib,iparm)=Q0(j,i,ib,1) 
227              enddo
228            enddo
229            replica(iparm)=replica(1)
230            umbrella(iparm)=umbrella(1)
231            read_iset(iparm)=read_iset(1)
232          enddo
233        enddo
234         
235       endif
236
237       return
238       end
239 c-----------------------------------------------------------------------------
240       subroutine read_protein_data(*)
241       implicit none
242       include "DIMENSIONS"
243       include "DIMENSIONS.ZSCOPT"
244       include "DIMENSIONS.FREE"
245 #ifdef MPI
246       include "mpif.h"
247       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
248       include "COMMON.MPI"
249 #endif
250       include "COMMON.CHAIN"
251       include "COMMON.IOUNITS"
252       include "COMMON.PROT"
253       include "COMMON.PROTFILES"
254       include "COMMON.NAMES"
255       include "COMMON.FREE"
256       include "COMMON.OBCINKA"
257       character*64 nazwa
258       character*16000 controlcard
259       integer i,ii,ib,iR,iparm,ilen,iroof,nthr,npars
260       external ilen,iroof
261       if (hamil_rep) then
262         npars=1
263       else
264         npars=nparmset
265       endif
266
267       do iparm=1,npars
268
269 C Read names of files with conformation data.
270       if (replica(iparm)) then
271         nthr = 1
272       else
273         nthr = nT_h(iparm)
274       endif
275       do ib=1,nthr
276       do ii=1,nRR(ib,iparm)
277       write (iout,*) "Parameter set",iparm," temperature",ib,
278      & " window",ii
279       call flush(iout)
280       call card_concat(controlcard,.true.) 
281       write (iout,*) controlcard(:ilen(controlcard))
282       call readi(controlcard,"NFILE_BIN",nfile_bin(ii,ib,iparm),0)
283       call readi(controlcard,"NFILE_ASC",nfile_asc(ii,ib,iparm),0)
284       call readi(controlcard,"NFILE_CX",nfile_cx(ii,ib,iparm),0)
285       call readi(controlcard,"REC_START",rec_start(ii,ib,iparm),1)
286       call readi(controlcard,"REC_END",rec_end(ii,ib,iparm),
287      & maxstr*isampl(iparm)+rec_start(ii,ib,iparm)-1)
288       call reada(controlcard,"TIME_START",
289      &  time_start_collect(ii,ib,iparm),0.0d0)
290       call reada(controlcard,"TIME_END",time_end_collect(ii,ib,iparm),
291      &  1.0d10)
292       write (iout,*) "rec_start",rec_start(ii,ib,iparm),
293      & " rec_end",rec_end(ii,ib,iparm)
294       write (iout,*) "time_start",time_start_collect(ii,ib,iparm),
295      & " time_end",time_end_collect(ii,ib,iparm)
296       call flush(iout)
297       if (replica(iparm)) then
298         call readi(controlcard,"TOTRAJ",totraj(ii,iparm),1)
299         write (iout,*) "Number of trajectories",totraj(ii,iparm)
300         call flush(iout)
301       endif
302       if (nfile_bin(ii,ib,iparm).lt.2 
303      &    .and. nfile_asc(ii,ib,iparm).eq.0
304      &    .and. nfile_cx(ii,ib,iparm).eq.0) then
305         write (iout,*) "Error - no action specified!"
306         return1
307       endif
308       if (nfile_bin(ii,ib,iparm).gt.0) then
309         call card_concat(controlcard,.false.)
310         call split_string(controlcard,protfiles(1,1,ii,ib,iparm),
311      &   maxfile_prot,nfile_bin(ii,ib,iparm))
312 #ifdef DEBUG
313         write(iout,*)"nfile_bin",nfile_bin(ii,ib,iparm)
314         write(iout,*) (protfiles(i,1,ii,ib,iparm),
315      &    i=1,nfile_bin(ii,ib,iparm))
316 #endif
317       endif
318       if (nfile_asc(ii,ib,iparm).gt.0) then
319         call card_concat(controlcard,.false.)
320         call split_string(controlcard,protfiles(1,2,ii,ib,iparm),
321      &   maxfile_prot,nfile_asc(ii,ib,iparm))
322 #ifdef DEBUG
323         write(iout,*) "nfile_asc(ii,ib,iparm)",nfile_asc(ii,ib,iparm)
324         write(iout,*) (protfiles(i,2,ii,ib,iparm),
325      &    i=1,nfile_asc(ii,ib,iparm))
326 #endif
327       else if (nfile_cx(ii,ib,iparm).gt.0) then
328         call card_concat(controlcard,.false.)
329         call split_string(controlcard,protfiles(1,2,ii,ib,iparm),
330      &   maxfile_prot,nfile_cx(ii,ib,iparm))
331 #ifdef DEBUG
332         write(iout,*) "nfile_cx(ii,ib,iparm)",nfile_cx(ii,ib,iparm)
333         write(iout,*) (protfiles(i,2,ii,ib,iparm),
334      &    i=1,nfile_cx(ii,ib,iparm))
335 #endif
336       endif
337       call flush(iout)
338       enddo
339       enddo
340
341       enddo
342
343       return
344       end
345 c-------------------------------------------------------------------------------
346       subroutine opentmp(islice,iunit,bprotfile_temp)
347       implicit none
348       include "DIMENSIONS"
349       include "DIMENSIONS.ZSCOPT"
350       include "DIMENSIONS.FREE"
351 #ifdef MPI
352       include "mpif.h"
353       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
354       include "COMMON.MPI"
355 #endif
356       include "COMMON.IOUNITS"
357       include "COMMON.PROTFILES"
358       include "COMMON.PROT"
359       include "COMMON.FREE"
360       character*64 bprotfile_temp
361       character*3 liczba,liczba2
362       character*2 liczba1
363       integer iunit,islice
364       integer ilen,iroof
365       external ilen,iroof
366       logical lerr
367
368       write (liczba1,'(bz,i2.2)') islice
369       write (liczba,'(bz,i3.3)') me
370 #ifdef MPI
371 c      write (iout,*) "separate_parset ",separate_parset,
372 c     &  " myparm",myparm
373       if (separate_parset) then
374       write (liczba2,'(bz,i3.3)') myparm
375       bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
376      &  prefix(:ilen(prefix))//liczba//"_"//liczba2//".xbin.tmp"//liczba1
377       open (iunit,file=bprotfile_temp,status="unknown",
378      &    form="unformatted",access="direct",recl=lenrec)
379       else
380       bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
381      &  prefix(:ilen(prefix))//liczba//".xbin.tmp"//liczba1
382       open (iunit,file=bprotfile_temp,status="unknown",
383      &    form="unformatted",access="direct",recl=lenrec)
384       endif
385 #else
386       bprotfile_temp = scratchdir(:ilen(scratchdir))//
387      &  "/"//prefix(:ilen(prefix))//".xbin.tmp"//liczba1
388       open (iunit,file=bprotfile_temp,status="unknown",
389      &    form="unformatted",access="direct",recl=lenrec)
390 #endif      
391 c      write (iout,*) "OpenTmp iunit",iunit," bprotfile_temp",
392 c     &  bprotfile_temp
393 c      call flush(iout)
394       return
395       end
396 c-------------------------------------------------------------------------------
397       subroutine read_database(*)
398       implicit none
399       include "DIMENSIONS"
400       include "DIMENSIONS.ZSCOPT"
401       include "DIMENSIONS.FREE"
402 #ifdef MPI
403       include "mpif.h"
404       integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
405       include "COMMON.MPI"
406 #endif
407       include "COMMON.CHAIN"
408       include "COMMON.IOUNITS"
409       include "COMMON.PROTFILES"
410       include "COMMON.NAMES"
411       include "COMMON.VAR"
412       include "COMMON.GEO"
413       include "COMMON.ENEPS"
414       include "COMMON.PROT"
415       include "COMMON.INTERACT"
416       include "COMMON.FREE"
417       include "COMMON.SBRIDGE"
418       include "COMMON.OBCINKA"
419       real*4 csingle(3,maxres2)
420       character*64 nazwa,bprotfile_temp
421       character*3 liczba
422       character*2 liczba1
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
427       integer ilen,iroof
428       external ilen,iroof
429       double precision rmsdev,energia(0:max_ene),efree,eini,temp
430       double precision prop(maxQ)
431       integer ntot_all(maxslice,0:maxprocs-1), maxslice_buff
432       integer iparm,ib,iib,ir,nprop,nthr,npars
433       double precision etot,time
434       integer ixdrf,iret 
435       logical lerr,linit
436
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
439       lenrec=lenrec2+8
440       write (iout,*) "lenrec",lenrec," lenrec1",lenrec1,
441      &  " lenrec2",lenrec2
442
443       do i=1,nQ
444         prop(i)=0.0d0
445       enddo
446       do islice=1,nslice
447         ll(islice)=0
448         mm(islice)=0
449       enddo
450       write (iout,*) "nparmset",nparmset
451       if (hamil_rep) then
452         npars=1
453       else
454         npars=nparmset
455       endif
456       do iparm=1,npars
457
458       if (replica(iparm)) then
459         nthr = 1
460       else
461         nthr = nT_h(iparm)
462       endif
463
464       do ib=1,nthr
465       do iR=1,nRR(ib,iparm)
466
467       write (iout,*) "iparm",iparm," ib",ib," iR",iR," nQ",nQ
468       do islice=1,nslice
469         jj(islice)=0
470         kk(islice)=0
471       enddo
472
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)
482           ii=0
483           do islice=1,nslice
484             call opentmp(islice,ientout,bprotfile_temp)
485             call bxread(nazwa,ii,jj(islice),kk(islice),ll(islice),
486      &        mm(islice),iR,ib,iparm)
487             close(ientout)
488           enddo
489           close(ientin)
490         enddo
491       ENDIF ! NFILE_BIN>0
492 c
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
495 c DA scratchfile.
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))
501           ii=0
502           call xread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm)
503         enddo ! if
504       ENDIF
505       IF (NFILE_CX(iR,ib,iparm).gt.0) THEN
506 c Read conformations from cx files and write them to a binary
507 c DA scratchfile.
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))
512           ii=0
513           print *,"Calling cxread"
514           call cxread(nazwa,ii,jj,kk,ll,mm,iR,ib,iparm,
515      &       *1111)
516           close(ientout)
517           write (iout,*) "exit cxread"
518           call flush(iout)
519         enddo
520       ENDIF
521
522       do islice=1,nslice
523         stot(islice)=stot(islice)+jj(islice)
524       enddo
525
526       enddo
527       enddo
528       write (iout,*) "IPARM",iparm
529       enddo
530
531       if (nslice.eq.1) then
532 #ifdef MPI
533         write (liczba,'(bz,i3.3)') me
534         bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
535      &    prefix(:ilen(prefix))//liczba//".xbin.tmp"
536 #else
537         bprotfile_temp = scratchdir(:ilen(scratchdir))//
538      &     "/"//prefix(:ilen(prefix))//".xbin.tmp"
539 #endif
540         write(iout,*) mm(1)," conformations read",ll(1),
541      &    " conformations written to ", 
542      &    bprotfile_temp(:ilen(bprotfile_temp))
543       else
544         do islice=1,nslice
545           write (liczba1,'(bz,i2.2)') islice
546 #ifdef MPI
547           write (liczba,'(bz,i3.3)') me
548           bprotfile_temp = scratchdir(:ilen(scratchdir))//"/"//
549      &      prefix(:ilen(prefix))//liczba//".xbin.tmp"//liczba1
550 #else
551           bprotfile_temp = scratchdir(:ilen(scratchdir))//
552      &       "/"//prefix(:ilen(prefix))//".xbin.tmp"//liczba1
553 #endif
554           write(iout,*) mm(islice)," conformations read",ll(islice),
555      &    " conformations written to ", 
556      &    bprotfile_temp(:ilen(bprotfile_temp))
557         enddo
558       endif
559
560 #ifdef MPI
561 c Check if everyone has the same number of conformations
562
563 c      call MPI_ALLgather(MPI_IN_PLACE,stot(1),MPI_DATATYPE_NULL,
564 c     &  ntot_all(1,0),maxslice,MPI_INTEGER,MPI_Comm_World,IERROR)
565
566       maxslice_buff=maxslice
567
568       call MPI_Allgather(stot(1),maxslice_buff,MPI_INTEGER,
569      &  ntot_all(1,0),maxslice,MPI_INTEGER,MPI_Comm_World,IERROR)
570       lerr=.false.
571       do i=0,nprocs-1
572         if (i.ne.me) then
573           do islice=1,nslice
574           if (stot(islice).ne.ntot_all(islice,i)) then
575             write (iout,*) "Number of conformations at processor",i,
576      &       " differs from that at processor",me,
577      &       stot(islice),ntot_all(islice,i)," slice",islice
578             lerr = .true.
579           endif
580           enddo
581         endif
582       enddo 
583       if (lerr) then
584         write (iout,*)
585         write (iout,*) "Numbers of conformations read by processors"
586         write (iout,*)
587         do i=0,nprocs-1
588           write (iout,'(8i10)') i,(ntot_all(islice,i),islice=1,nslice)
589         enddo
590         write (iout,*) "Calculation terminated."
591         call flush(iout)
592         return1
593       endif
594       do islice=1,nslice
595         ntot(islice)=stot(islice)
596       enddo
597       return
598 #endif
599  1111 write(iout,*) "Error opening coordinate file ",nazwa(:ilen(nazwa))
600       call flush(iout)
601       return1
602       end
603 c------------------------------------------------------------------------------
604       subroutine card_concat(card,to_upper)
605       implicit none
606       include 'DIMENSIONS.ZSCOPT'
607       include "COMMON.IOUNITS"
608       character*(*) card
609       character*80 karta,ucase
610       logical to_upper
611       integer ilen
612       external ilen
613       read (inp,'(a)') karta
614       if (to_upper) karta=ucase(karta)
615       card=' '
616       do while (karta(80:80).eq.'&')
617         card=card(:ilen(card)+1)//karta(:79)
618         read (inp,'(a)') karta
619         if (to_upper) karta=ucase(karta)
620       enddo
621       card=card(:ilen(card)+1)//karta
622       return
623       end
624 c------------------------------------------------------------------------------
625       subroutine readi(rekord,lancuch,wartosc,default)
626       implicit none
627       character*(*) rekord,lancuch
628       integer wartosc,default
629       integer ilen,iread
630       external ilen
631       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
632       if (iread.eq.0) then
633         wartosc=default
634         return
635       endif
636       iread=iread+ilen(lancuch)+1
637       read (rekord(iread:),*) wartosc
638       return
639       end
640 c----------------------------------------------------------------------------
641       subroutine reada(rekord,lancuch,wartosc,default)
642       implicit none
643       character*(*) rekord,lancuch
644       character*80 aux
645       double precision wartosc,default
646       integer ilen,iread
647       external ilen
648       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
649       if (iread.eq.0) then
650         wartosc=default
651         return
652       endif
653       iread=iread+ilen(lancuch)+1
654       read (rekord(iread:),*) wartosc
655       return
656       end
657 c----------------------------------------------------------------------------
658       subroutine multreadi(rekord,lancuch,tablica,dim,default)
659       implicit none
660       integer dim,i
661       integer tablica(dim),default
662       character*(*) rekord,lancuch
663       character*80 aux
664       integer ilen,iread
665       external ilen
666       do i=1,dim
667         tablica(i)=default
668       enddo
669       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
670       if (iread.eq.0) return
671       iread=iread+ilen(lancuch)+1
672       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
673    10 return
674       end
675 c----------------------------------------------------------------------------
676       subroutine multreada(rekord,lancuch,tablica,dim,default)
677       implicit none
678       integer dim,i
679       double precision tablica(dim),default
680       character*(*) rekord,lancuch
681       character*80 aux
682       integer ilen,iread
683       external ilen
684       do i=1,dim
685         tablica(i)=default
686       enddo
687       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
688       if (iread.eq.0) return
689       iread=iread+ilen(lancuch)+1
690       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
691    10 return
692       end
693 c----------------------------------------------------------------------------
694       subroutine reads(rekord,lancuch,wartosc,default)
695       implicit none
696       character*(*) rekord,lancuch,wartosc,default
697       character*80 aux
698       integer ilen,lenlan,lenrec,iread,ireade
699       external ilen
700       logical iblnk
701       external iblnk
702       lenlan=ilen(lancuch)
703       lenrec=ilen(rekord)
704       iread=index(rekord,lancuch(:lenlan)//"=")
705 c      print *,"rekord",rekord," lancuch",lancuch
706 c      print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
707       if (iread.eq.0) then
708         wartosc=default
709         return
710       endif
711       iread=iread+lenlan+1
712 c      print *,"iread",iread
713 c      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
714       do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
715         iread=iread+1
716 c      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
717       enddo
718 c      print *,"iread",iread
719       if (iread.gt.lenrec) then
720          wartosc=default
721         return
722       endif
723       ireade=iread+1
724 c      print *,"ireade",ireade
725       do while (ireade.lt.lenrec .and.
726      &   .not.iblnk(rekord(ireade:ireade)))
727         ireade=ireade+1
728       enddo
729       wartosc=rekord(iread:ireade)
730       return
731       end
732 c----------------------------------------------------------------------------
733       subroutine multreads(rekord,lancuch,tablica,dim,default)
734       implicit none
735       integer dim,i
736       character*(*) rekord,lancuch,tablica(dim),default
737       character*80 aux
738       integer ilen,lenlan,lenrec,iread,ireade
739       external ilen
740       logical iblnk
741       external iblnk
742       do i=1,dim
743         tablica(i)=default
744       enddo
745       lenlan=ilen(lancuch)
746       lenrec=ilen(rekord)
747       iread=index(rekord,lancuch(:lenlan)//"=")
748 c      print *,"rekord",rekord," lancuch",lancuch
749 c      print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
750       if (iread.eq.0) return
751       iread=iread+lenlan+1
752       do i=1,dim
753 c      print *,"iread",iread
754 c      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
755       do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
756         iread=iread+1
757 c      print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
758       enddo
759 c      print *,"iread",iread
760       if (iread.gt.lenrec) return
761       ireade=iread+1
762 c      print *,"ireade",ireade
763       do while (ireade.lt.lenrec .and.
764      &   .not.iblnk(rekord(ireade:ireade)))
765         ireade=ireade+1
766       enddo
767       tablica(i)=rekord(iread:ireade)
768       iread=ireade+1
769       enddo
770       end
771 c----------------------------------------------------------------------------
772       subroutine split_string(rekord,tablica,dim,nsub)
773       implicit none
774       integer dim,nsub,i,ii,ll,kk
775       character*(*) tablica(dim)
776       character*(*) rekord
777       integer ilen
778       external ilen
779       do i=1,dim
780         tablica(i)=" "
781       enddo
782       ii=1
783       ll = ilen(rekord)
784       nsub=0
785       do i=1,dim
786 C Find the start of term name
787         kk = 0
788         do while (ii.le.ll .and. rekord(ii:ii).eq." ") 
789           ii = ii+1
790         enddo
791 C Parse the name into TABLICA(i) until blank found
792         do while (ii.le.ll .and. rekord(ii:ii).ne." ") 
793           kk = kk+1
794           tablica(i)(kk:kk)=rekord(ii:ii)
795           ii = ii+1
796         enddo
797         if (kk.gt.0) nsub=nsub+1
798         if (ii.gt.ll) return
799       enddo
800       return
801       end
802 c--------------------------------------------------------------------------------
803       integer function iroof(n,m)
804       ii = n/m
805       if (ii*m .lt. n) ii=ii+1
806       iroof = ii
807       return
808       end