D i L aminokwasy dzialajace dla src_CSA
[unres.git] / source / unres / src_CSA_DiL / readrtns_csa.F
1       subroutine readrtns
2       implicit real*8 (a-h,o-z)
3       include 'DIMENSIONS'
4 #ifdef MPI
5       include 'mpif.h'
6 #endif
7       include 'COMMON.SETUP'
8       include 'COMMON.CONTROL'
9       include 'COMMON.SBRIDGE'
10       include 'COMMON.IOUNITS'
11       logical file_exist
12 C Read force-field parameters except weights
13       call parmread
14 C Read job setup parameters
15       call read_control
16 C Read control parameters for energy minimzation if required
17       if (minim) call read_minim
18 C Read MCM control parameters if required
19 c      if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
20 C Read MD control parameters if reqjuired
21 c      if (modecalc.eq.12) call read_MDpar
22 C Read MREMD control parameters if required
23 c      if (modecalc.eq.14) then 
24 c         call read_MDpar
25 c         call read_REMDpar
26 c      endif
27 C Read MUCA control parameters if required
28 c      if (lmuca) call read_muca
29 C Read CSA control parameters if required (from fort.40 if exists
30 C otherwise from general input file)
31       if (modecalc.eq.8) then
32        inquire (file="fort.40",exist=file_exist)
33        if (.not.file_exist) call csaread
34       endif 
35 cfmc      if (modecalc.eq.10) call mcmfread
36 C Read molecule information, molecule geometry, energy-term weights, and
37 C restraints if requested
38       call molread
39       
40 C Print restraint information
41 #ifdef MPI
42       if (.not. out1file .or. me.eq.king) then
43 #endif
44       if (nhpb.gt.nss) 
45      &write (iout,'(a,i5,a)') "The following",nhpb-nss,
46      & " distance constraints have been imposed"
47       do i=nss+1,nhpb
48         write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
49       enddo
50 #ifdef MPI
51       endif
52 #endif
53 c      print *,"Processor",myrank," leaves READRTNS"
54       return
55       end
56 C-------------------------------------------------------------------------------
57       subroutine read_control
58 C
59 C Read contorl data
60 C
61       implicit real*8 (a-h,o-z)
62       include 'DIMENSIONS'
63 #ifdef MP
64       include 'mpif.h'
65       logical OKRandom, prng_restart
66       real*8  r1
67 #endif
68       include 'COMMON.IOUNITS'
69       include 'COMMON.TIME1'
70 c      include 'COMMON.THREAD'
71       include 'COMMON.SBRIDGE'
72       include 'COMMON.CONTROL'
73       include 'COMMON.MCM'
74 c      include 'COMMON.MAP'
75       include 'COMMON.HEADER'
76 c      include 'COMMON.CSA'
77       include 'COMMON.CHAIN'
78 c      include 'COMMON.MUCA'
79 c      include 'COMMON.MD'
80       include 'COMMON.FFIELD'
81       include 'COMMON.SETUP'
82       COMMON /MACHSW/ KDIAG,ICORFL,IXDR
83       character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
84       character*80 ucase
85       character*320 controlcard
86
87       nglob_csa=0
88       eglob_csa=1d99
89       nmin_csa=0
90       read (INP,'(a)') titel
91       call card_concat(controlcard)
92 c      out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
93 c      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
94       call reada(controlcard,'SEED',seed,0.0D0)
95       call random_init(seed)
96 C Set up the time limit (caution! The time must be input in minutes!)
97       read_cart=index(controlcard,'READ_CART').gt.0
98       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
99       call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
100       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
101       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
102       call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
103       call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
104       call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
105       call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
106       call reada(controlcard,'DRMS',drms,0.1D0)
107       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
108        write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc 
109        write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1 
110        write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max 
111        write (iout,'(a,f10.1)')'DRMS    = ',drms 
112        write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm 
113        write (iout,'(a,f10.1)') 'Time limit (min):',timlim
114       endif
115       call readi(controlcard,'NZ_START',nz_start,0)
116       call readi(controlcard,'NZ_END',nz_end,0)
117       call readi(controlcard,'IZ_SC',iz_sc,0)
118       timlim=60.0D0*timlim
119       safety = 60.0d0*safety
120       timem=timlim
121       modecalc=0
122       call reada(controlcard,"T_BATH",t_bath,300.0d0)
123       minim=(index(controlcard,'MINIMIZE').gt.0)
124       dccart=(index(controlcard,'CART').gt.0)
125       overlapsc=(index(controlcard,'OVERLAP').gt.0)
126       overlapsc=.not.overlapsc
127       searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
128       searchsc=.not.searchsc
129       sideadd=(index(controlcard,'SIDEADD').gt.0)
130       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
131       outpdb=(index(controlcard,'PDBOUT').gt.0)
132       outmol2=(index(controlcard,'MOL2OUT').gt.0)
133       pdbref=(index(controlcard,'PDBREF').gt.0)
134       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
135       indpdb=index(controlcard,'PDBSTART')
136       extconf=(index(controlcard,'EXTCONF').gt.0)
137       call readi(controlcard,'IPRINT',iprint,0)
138       call readi(controlcard,'MAXGEN',maxgen,10000)
139       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
140       call readi(controlcard,"KDIAG",kdiag,0)
141       call readi(controlcard,"RESCALE_MODE",rescale_mode,0)
142       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
143      & write (iout,*) "RESCALE_MODE",rescale_mode
144       split_ene=index(controlcard,'SPLIT_ENE').gt.0
145       if (index(controlcard,'REGULAR').gt.0.0D0) then
146         call reada(controlcard,'WEIDIS',weidis,0.1D0)
147         modecalc=1
148         refstr=.true.
149       endif
150       if (index(controlcard,'CHECKGRAD').gt.0) then
151         modecalc=5
152         if (index(controlcard,'CART').gt.0) then
153           icheckgrad=1
154         elseif (index(controlcard,'CARINT').gt.0) then
155           icheckgrad=2
156         else
157           icheckgrad=3
158         endif
159       elseif (index(controlcard,'THREAD').gt.0) then
160         modecalc=2
161         call readi(controlcard,'THREAD',nthread,0)
162         if (nthread.gt.0) then
163           call reada(controlcard,'WEIDIS',weidis,0.1D0)
164         else
165           if (fg_rank.eq.0)
166      &    write (iout,'(a)')'A number has to follow the THREAD keyword.'
167           stop 'Error termination in Read_Control.'
168         endif
169       else if (index(controlcard,'MCMA').gt.0) then
170         modecalc=3
171       else if (index(controlcard,'MCEE').gt.0) then
172         modecalc=6
173       else if (index(controlcard,'MULTCONF').gt.0) then
174         modecalc=4
175       else if (index(controlcard,'MAP').gt.0) then
176         modecalc=7
177         call readi(controlcard,'MAP',nmap,0)
178       else if (index(controlcard,'CSA').gt.0) then
179         modecalc=8
180 crc      else if (index(controlcard,'ZSCORE').gt.0) then
181 crc   
182 crc  ZSCORE is rm from UNRES, modecalc=9 is available
183 crc
184 crc        modecalc=9
185 cfcm      else if (index(controlcard,'MCMF').gt.0) then
186 cfmc        modecalc=10
187       else if (index(controlcard,'SOFTREG').gt.0) then
188         modecalc=11
189       else if (index(controlcard,'CHECK_BOND').gt.0) then
190         modecalc=-1
191       else if (index(controlcard,'TEST').gt.0) then
192         modecalc=-2
193       else if (index(controlcard,'MD').gt.0) then
194         modecalc=12
195       else if (index(controlcard,'RE ').gt.0) then
196         modecalc=14
197       endif
198
199       lmuca=index(controlcard,'MUCA').gt.0
200       call readi(controlcard,'MUCADYN',mucadyn,0)      
201       call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
202       if (lmuca .and. (me.eq.king .or. .not.out1file )) 
203      & then
204        write (iout,*) 'MUCADYN=',mucadyn
205        write (iout,*) 'MUCASMOOTH=',muca_smooth
206       endif
207
208       iscode=index(controlcard,'ONE_LETTER')
209       indphi=index(controlcard,'PHI')
210       indback=index(controlcard,'BACK')
211       iranconf=index(controlcard,'RAND_CONF')
212       i2ndstr=index(controlcard,'USE_SEC_PRED')
213       gradout=index(controlcard,'GRADOUT').gt.0
214       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
215       
216       if(me.eq.king.or..not.out1file)
217      & write (iout,'(2a)') diagmeth(kdiag),
218      &  ' routine used to diagonalize matrices.'
219       return
220       end
221 c------------------------------------------------------------------------------
222       subroutine molread
223 C
224 C Read molecular data.
225 C
226       implicit real*8 (a-h,o-z)
227       include 'DIMENSIONS'
228 #ifdef MPI
229       include 'mpif.h'
230       integer error_msg
231 #endif
232       include 'COMMON.IOUNITS'
233       include 'COMMON.GEO'
234       include 'COMMON.VAR'
235       include 'COMMON.INTERACT'
236       include 'COMMON.LOCAL'
237       include 'COMMON.NAMES'
238       include 'COMMON.CHAIN'
239       include 'COMMON.FFIELD'
240       include 'COMMON.SBRIDGE'
241       include 'COMMON.HEADER'
242       include 'COMMON.CONTROL'
243 c      include 'COMMON.DBASE'
244 c      include 'COMMON.THREAD'
245       include 'COMMON.CONTACTS'
246       include 'COMMON.TORCNSTR'
247       include 'COMMON.TIME1'
248       include 'COMMON.BOUNDS'
249 c      include 'COMMON.MD'
250 c      include 'COMMON.REMD'
251       include 'COMMON.SETUP'
252       character*4 sequence(maxres)
253       integer rescode
254       double precision x(maxvar)
255       character*256 pdbfile
256       character*320 weightcard
257       character*80 weightcard_t,ucase
258       dimension itype_pdb(maxres)
259       common /pizda/ itype_pdb
260       logical seq_comp,fail
261       double precision energia(0:n_ene)
262       integer ilen
263       external ilen
264 C
265 C Body
266 C
267 C Read weights of the subsequent energy terms.
268
269        call card_concat(weightcard)
270        call reada(weightcard,'WLONG',wlong,1.0D0)
271        call reada(weightcard,'WSC',wsc,wlong)
272        call reada(weightcard,'WSCP',wscp,wlong)
273        call reada(weightcard,'WELEC',welec,1.0D0)
274        call reada(weightcard,'WVDWPP',wvdwpp,welec)
275        call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
276        call reada(weightcard,'WCORR4',wcorr4,0.0D0)
277        call reada(weightcard,'WCORR5',wcorr5,0.0D0)
278        call reada(weightcard,'WCORR6',wcorr6,0.0D0)
279        call reada(weightcard,'WTURN3',wturn3,1.0D0)
280        call reada(weightcard,'WTURN4',wturn4,1.0D0)
281        call reada(weightcard,'WTURN6',wturn6,1.0D0)
282        call reada(weightcard,'WSCCOR',wsccor,1.0D0)
283        call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
284        call reada(weightcard,'WBOND',wbond,1.0D0)
285        call reada(weightcard,'WTOR',wtor,1.0D0)
286        call reada(weightcard,'WTORD',wtor_d,1.0D0)
287        call reada(weightcard,'WANG',wang,1.0D0)
288        call reada(weightcard,'WSCLOC',wscloc,1.0D0)
289 C     Juyong
290        call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
291        call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
292        call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
293        call reada(weightcard,'WDFAB',wdfa_beta,0.0d0)
294 C       
295        call reada(weightcard,'SCAL14',scal14,0.4D0)
296        call reada(weightcard,'SCALSCP',scalscp,1.0d0)
297        call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
298        call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
299        call reada(weightcard,'TEMP0',temp0,300.0d0)
300        if (index(weightcard,'SOFT').gt.0) ipot=6
301 C 12/1/95 Added weight for the multi-body term WCORR
302        call reada(weightcard,'WCORRH',wcorr,1.0D0)
303        if (wcorr4.gt.0.0d0) wcorr=wcorr4
304        weights(1)=wsc
305        weights(2)=wscp
306        weights(3)=welec
307        weights(4)=wcorr
308        weights(5)=wcorr5
309        weights(6)=wcorr6
310        weights(7)=wel_loc
311        weights(8)=wturn3
312        weights(9)=wturn4
313        weights(10)=wturn6
314        weights(11)=wang
315        weights(12)=wscloc
316        weights(13)=wtor
317        weights(14)=wtor_d
318        weights(15)=wstrain
319        weights(16)=wvdwpp
320        weights(17)=wbond
321        weights(18)=scal14
322        weights(21)=wsccor
323 C     JUYONG
324        weights(24)=wdfa_dist
325        weights(25)=wdfa_tor
326        weights(26)=wdfa_nei
327        weights(27)=wdfa_beta
328 C
329
330       if(me.eq.king.or..not.out1file)
331      & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
332      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
333      &  wturn4,wturn6,
334      &  wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
335
336    10 format (/'Energy-term weights (unscaled):'//
337      & 'WSCC=   ',f10.6,' (SC-SC)'/
338      & 'WSCP=   ',f10.6,' (SC-p)'/
339      & 'WELEC=  ',f10.6,' (p-p electr)'/
340      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
341      & 'WBOND=  ',f10.6,' (stretching)'/
342      & 'WANG=   ',f10.6,' (bending)'/
343      & 'WSCLOC= ',f10.6,' (SC local)'/
344      & 'WTOR=   ',f10.6,' (torsional)'/
345      & 'WTORD=  ',f10.6,' (double torsional)'/
346      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
347      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
348      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
349      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
350      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
351      & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
352      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
353      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
354      & 'WTURN6= ',f10.6,' (turns, 6th order)'/
355      & 'WDFA_D= ',f10.6,' (DFA, distance)'   /
356      & 'WDFA_T= ',f10.6,' (DFA, torsional)'   /
357      & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)'   /
358      & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
359
360       if(me.eq.king.or..not.out1file)then
361        if (wcorr4.gt.0.0d0) then
362         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
363      &   'between contact pairs of peptide groups'
364         write (iout,'(2(a,f5.3/))') 
365      &  'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
366      &  'Range of quenching the correlation terms:',2*delt_corr 
367        else if (wcorr.gt.0.0d0) then
368         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
369      &   'between contact pairs of peptide groups'
370        endif
371        write (iout,'(a,f8.3)') 
372      &  'Scaling factor of 1,4 SC-p interactions:',scal14
373        write (iout,'(a,f8.3)') 
374      &  'General scaling factor of SC-p interactions:',scalscp
375       endif
376       r0_corr=cutoff_corr-delt_corr
377       do i=1,20
378         aad(i,1)=scalscp*aad(i,1)
379         aad(i,2)=scalscp*aad(i,2)
380         bad(i,1)=scalscp*bad(i,1)
381         bad(i,2)=scalscp*bad(i,2)
382       enddo
383 c      call rescale_weights(t_bath)
384       if(me.eq.king.or..not.out1file)
385      & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
386      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
387      &  wturn4,wturn6,
388      &  wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
389
390    22 format (/'Energy-term weights (scaled):'//
391      & 'WSCC=   ',f10.6,' (SC-SC)'/
392      & 'WSCP=   ',f10.6,' (SC-p)'/
393      & 'WELEC=  ',f10.6,' (p-p electr)'/
394      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
395      & 'WBOND=  ',f10.6,' (stretching)'/
396      & 'WANG=   ',f10.6,' (bending)'/
397      & 'WSCLOC= ',f10.6,' (SC local)'/
398      & 'WTOR=   ',f10.6,' (torsional)'/
399      & 'WTORD=  ',f10.6,' (double torsional)'/
400      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
401      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
402      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
403      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
404      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
405      & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
406      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
407      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
408      & 'WTURN6= ',f10.6,' (turns, 6th order)'/
409      & 'WDFA_D= ',f10.6,' (DFA, distance)'   /
410      & 'WDFA_T= ',f10.6,' (DFA, torsional)'   /
411      & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)'   /
412      & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
413
414       if(me.eq.king.or..not.out1file)
415      & write (iout,*) "Reference temperature for weights calculation:",
416      &  temp0
417       call reada(weightcard,"D0CM",d0cm,3.78d0)
418       call reada(weightcard,"AKCM",akcm,15.1d0)
419       call reada(weightcard,"AKTH",akth,11.0d0)
420       call reada(weightcard,"AKCT",akct,12.0d0)
421       call reada(weightcard,"V1SS",v1ss,-1.08d0)
422       call reada(weightcard,"V2SS",v2ss,7.61d0)
423       call reada(weightcard,"V3SS",v3ss,13.7d0)
424       call reada(weightcard,"EBR",ebr,-5.50D0)
425       if(me.eq.king.or..not.out1file) then
426        write (iout,*) "Parameters of the SS-bond potential:"
427        write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
428      & " AKCT",akct
429        write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
430        write (iout,*) "EBR",ebr
431        print *,'indpdb=',indpdb,' pdbref=',pdbref
432       endif
433       if (indpdb.gt.0 .or. pdbref) then
434         read(inp,'(a)') pdbfile
435         if(me.eq.king.or..not.out1file)
436      &   write (iout,'(2a)') 'PDB data will be read from file ',
437      &   pdbfile(:ilen(pdbfile))
438         open(ipdbin,file=pdbfile,status='old',err=33)
439         goto 34 
440   33    write (iout,'(a)') 'Error opening PDB file.'
441         stop
442   34    continue
443 c        print *,'Begin reading pdb data'
444         call readpdb
445 c        print *,'Finished reading pdb data'
446         if(me.eq.king.or..not.out1file)
447      &   write (iout,'(a,i3,a,i3)')'nsup=',nsup,
448      &   ' nstart_sup=',nstart_sup
449         do i=1,nres
450           itype_pdb(i)=itype(i)
451         enddo
452         close (ipdbin)
453         nnt=nstart_sup
454         nct=nstart_sup+nsup-1
455         call contact(.false.,ncont_ref,icont_ref,co)
456
457         if (sideadd) then 
458          if(me.eq.king.or..not.out1file)
459      &    write(iout,*)'Adding sidechains'
460          maxsi=1000
461          do i=2,nres-1
462           iti=itype(i)
463           if (iti.ne.10) then
464             nsi=0
465             fail=.true.
466             do while (fail.and.nsi.le.maxsi)
467 c              call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
468               nsi=nsi+1
469             enddo
470             if(fail) write(iout,*)'Adding sidechain failed for res ',
471      &              i,' after ',nsi,' trials'
472           endif
473          enddo
474         endif  
475       endif
476
477       if (indpdb.eq.0) then
478 C Read sequence if not taken from the pdb file.
479         read (inp,*) nres
480 c        print *,'nres=',nres
481         if (iscode.gt.0) then
482           read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
483         else
484           read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
485         endif
486 C Convert sequence to numeric code
487         do i=1,nres
488           itype(i)=rescode(i,sequence(i),iscode)
489         enddo
490 C Assign initial virtual bond lengths
491         do i=2,nres
492           vbld(i)=vbl
493           vbld_inv(i)=vblinv
494         enddo
495         do i=2,nres-1
496           vbld(i+nres)=dsc(itype(i))
497           vbld_inv(i+nres)=dsc_inv(itype(i))
498 c          write (iout,*) "i",i," itype",itype(i),
499 c     &      " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
500         enddo
501       endif 
502 c      print *,nres
503 c      print '(20i4)',(itype(i),i=1,nres)
504       do i=1,nres
505 #ifdef PROCOR
506         if (itype(i).eq.21 .or. itype(i+1).eq.21) then
507 #else
508         if (itype(i).eq.21) then
509 #endif
510           itel(i)=0
511 #ifdef PROCOR
512         else if (iabs(itype(i+1)).ne.20) then
513 #else
514         else if (iabs(itype(i)).ne.20) then
515 #endif
516           itel(i)=1
517         else
518           itel(i)=2
519         endif  
520       enddo
521       if(me.eq.king.or..not.out1file)then
522        write (iout,*) "ITEL"
523        do i=1,nres-1
524          write (iout,*) i,itype(i),itel(i)
525        enddo
526        print *,'Call Read_Bridge.'
527       endif
528       call read_bridge
529 C 8/13/98 Set limits to generating the dihedral angles
530       do i=1,nres
531         phibound(1,i)=-pi
532         phibound(2,i)=pi
533       enddo
534       read (inp,*) ndih_constr
535       if (ndih_constr.gt.0) then
536         read (inp,*) ftors
537         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
538         if(me.eq.king.or..not.out1file)then
539          write (iout,*) 
540      &   'There are',ndih_constr,' constraints on phi angles.'
541          do i=1,ndih_constr
542           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
543          enddo
544         endif
545         do i=1,ndih_constr
546           phi0(i)=deg2rad*phi0(i)
547           drange(i)=deg2rad*drange(i)
548         enddo
549         if(me.eq.king.or..not.out1file)
550      &   write (iout,*) 'FTORS',ftors
551         do i=1,ndih_constr
552           ii = idih_constr(i)
553           phibound(1,ii) = phi0(i)-drange(i)
554           phibound(2,ii) = phi0(i)+drange(i)
555         enddo 
556       endif
557       nnt=1
558 #ifdef MPI
559       if (me.eq.king) then
560 #endif
561        write (iout,'(a)') 'Boundaries in phi angle sampling:'
562        do i=1,nres
563          write (iout,'(a3,i5,2f10.1)') 
564      &   restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
565        enddo
566 #ifdef MP
567       endif
568 #endif
569       nct=nres
570 cd      print *,'NNT=',NNT,' NCT=',NCT
571       if (itype(1).eq.21) nnt=2
572       if (itype(nres).eq.21) nct=nct-1
573
574 C     Juyong:READ init_vars
575 C     Initialize variables!
576 C     Juyong:READ read_info
577 C     READ fragment information!!
578 C     both routines should be in dfa.F file!!
579
580       if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
581      &            wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
582        call init_dfa_vars
583        print*, 'init_dfa_vars finished!'
584        call read_dfa_info
585        print*, 'read_dfa_info finished!'
586       endif
587 C
588 C
589
590
591       if (pdbref) then
592         if(me.eq.king.or..not.out1file)
593      &   write (iout,'(a,i3)') 'nsup=',nsup
594         nstart_seq=nnt
595         if (nsup.le.(nct-nnt+1)) then
596           do i=0,nct-nnt+1-nsup
597             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
598               nstart_seq=nnt+i
599               goto 111
600             endif
601           enddo
602           write (iout,'(a)') 
603      &            'Error - sequences to be superposed do not match.'
604           stop
605         else
606           do i=0,nsup-(nct-nnt+1)
607             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1)) 
608      &      then
609               nstart_sup=nstart_sup+i
610               nsup=nct-nnt+1
611               goto 111
612             endif
613           enddo 
614           write (iout,'(a)') 
615      &            'Error - sequences to be superposed do not match.'
616         endif
617   111   continue
618         if (nsup.eq.0) nsup=nct-nnt
619         if (nstart_sup.eq.0) nstart_sup=nnt
620         if (nstart_seq.eq.0) nstart_seq=nnt
621         if(me.eq.king.or..not.out1file)  
622      &   write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
623      &                 ' nstart_seq=',nstart_seq
624       endif
625 c--- Zscore rms -------
626       if (nz_start.eq.0) nz_start=nnt
627       if (nz_end.eq.0 .and. nsup.gt.0) then
628         nz_end=nnt+nsup-1
629       else if (nz_end.eq.0) then
630         nz_end=nct
631       endif
632       if(me.eq.king.or..not.out1file)then
633        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
634        write (iout,*) 'IZ_SC=',iz_sc
635       endif
636 c----------------------
637       call init_int_table
638       if (refstr) then
639         if (.not.pdbref) then
640           call read_angles(inp,*38)
641           goto 39
642    38     write (iout,'(a)') 'Error reading reference structure.'
643 #ifdef MPI
644           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
645           stop 'Error reading reference structure'
646 #endif
647    39     call chainbuild
648           call setup_var
649 czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
650           nstart_sup=nnt
651           nstart_seq=nnt
652           nsup=nct-nnt+1
653           do i=1,2*nres
654             do j=1,3
655               cref(j,i)=c(j,i)
656             enddo
657           enddo
658           call contact(.true.,ncont_ref,icont_ref,co)
659         endif
660 c        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
661         call flush(iout)
662         if (constr_dist.gt.0) call read_dist_constr
663 c        write (iout,*) "After read_dist_constr nhpb",nhpb
664         call hpb_partition
665         if(me.eq.king.or..not.out1file)
666      &   write (iout,*) 'Contact order:',co
667         if (pdbref) then
668         if(me.eq.king.or..not.out1file)
669      &   write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
670         do i=1,ncont_ref
671           do j=1,2
672             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
673           enddo
674           if(me.eq.king.or..not.out1file)
675      &     write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
676      &     icont_ref(1,i),' ',
677      &     restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
678         enddo
679         endif
680       endif
681       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
682      &    .and. modecalc.ne.8 .and. modecalc.ne.9 .and. 
683      &    modecalc.ne.10) then
684 C If input structure hasn't been supplied from the PDB file read or generate
685 C initial geometry.
686         if (iranconf.eq.0 .and. .not. extconf) then
687           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
688      &     write (iout,'(a)') 'Initial geometry will be read in.'
689           if (read_cart) then
690             read(inp,'(8f10.5)',end=36,err=36)
691      &       ((c(l,k),l=1,3),k=1,nres),
692      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
693             call int_from_cart1(.false.)
694             do i=1,nres-1
695               do j=1,3
696                 dc(j,i)=c(j,i+1)-c(j,i)
697                 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
698               enddo
699             enddo
700             do i=nnt,nct
701               if (itype(i).ne.10) then
702                 do j=1,3
703                   dc(j,i+nres)=c(j,i+nres)-c(j,i) 
704                   dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
705                 enddo
706               endif
707             enddo
708             return
709           else
710             call read_angles(inp,*36)
711           endif
712           goto 37
713    36     write (iout,'(a)') 'Error reading angle file.'
714 #ifdef MPI
715           call mpi_finalize( MPI_COMM_WORLD,IERR )
716 #endif
717           stop 'Error reading angle file.'
718    37     continue 
719         else if (extconf) then
720          if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
721      &    write (iout,'(a)') 'Extended chain initial geometry.'
722          do i=3,nres
723           theta(i)=90d0*deg2rad
724          enddo
725          do i=4,nres
726           phi(i)=180d0*deg2rad
727          enddo
728          do i=2,nres-1
729           alph(i)=110d0*deg2rad
730          enddo
731          do i=2,nres-1
732           omeg(i)=-120d0*deg2rad
733          enddo
734         else
735           if(me.eq.king.or..not.out1file)
736      &     write (iout,'(a)') 'Random-generated initial geometry.'
737
738
739 #ifdef MPI
740           if (me.eq.king  .or. fg_rank.eq.0 .and. (
741      &           modecalc.eq.12 .or. modecalc.eq.14) ) then  
742 #endif
743             do itrial=1,100
744               itmp=1
745 c              call gen_rand_conf(itmp,*30)
746               goto 40
747    30         write (iout,*) 'Failed to generate random conformation',
748      &          ', itrial=',itrial
749               write (*,*) 'Processor:',me,
750      &          ' Failed to generate random conformation',
751      &          ' itrial=',itrial
752               call intout
753
754 #ifdef AIX
755               call flush_(iout)
756 #else
757               call flush(iout)
758 #endif
759             enddo
760             write (iout,'(a,i3,a)') 'Processor:',me,
761      &        ' error in generating random conformation.'
762             write (*,'(a,i3,a)') 'Processor:',me,
763      &        ' error in generating random conformation.'
764             call flush(iout)
765 #ifdef MPI
766             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)            
767    40       continue
768           endif
769 #else
770           do itrial=1,100
771             itmp=1
772 c            call gen_rand_conf(itmp,*31)
773             goto 40
774    31       write (iout,*) 'Failed to generate random conformation',
775      &        ', itrial=',itrial
776             write (*,*) 'Failed to generate random conformation',
777      &        ', itrial=',itrial
778           enddo
779           write (iout,'(a,i3,a)') 'Processor:',me,
780      &      ' error in generating random conformation.'
781           write (*,'(a,i3,a)') 'Processor:',me,
782      &      ' error in generating random conformation.'
783           stop
784    40     continue
785 #endif
786         endif
787       elseif (modecalc.eq.4) then
788         read (inp,'(a)') intinname
789         open (intin,file=intinname,status='old',err=333)
790         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
791      &  write (iout,'(a)') 'intinname',intinname
792         write (*,'(a)') 'Processor',myrank,' intinname',intinname
793         goto 334
794   333   write (iout,'(2a)') 'Error opening angle file ',intinname
795 #ifdef MPI 
796         call MPI_Finalize(MPI_COMM_WORLD,IERR)
797 #endif   
798         stop 'Error opening angle file.' 
799   334   continue
800
801       endif 
802 C Generate distance constraints, if the PDB structure is to be regularized. 
803       if (nthread.gt.0) then
804         call read_threadbase
805       endif
806       call setup_var
807       if (me.eq.king .or. .not. out1file)
808      & call intout
809       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
810         write (iout,'(/a,i3,a)') 
811      &  'The chain contains',ns,' disulfide-bridging cysteines.'
812         write (iout,'(20i4)') (iss(i),i=1,ns)
813         write (iout,'(/a/)') 'Pre-formed links are:' 
814         do i=1,nss
815           i1=ihpb(i)-nres
816           i2=jhpb(i)-nres
817           it1=itype(i1)
818           it2=itype(i2)
819           if (me.eq.king.or..not.out1file)
820      &    write (iout,'(2a,i3,3a,i3,a,3f10.3)')
821      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
822      &    ebr,forcon(i)
823         enddo
824         write (iout,'(a)')
825       endif
826 c       if (i2ndstr.gt.0) call secstrp2dihc
827 c      call geom_to_var(nvar,x)
828 c      call etotal(energia(0))
829 c      call enerprint(energia(0))
830 c      call briefout(0,etot)
831 c      stop
832 cd    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
833 cd    write (iout,'(a)') 'Variable list:'
834 cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
835 #ifdef MPI
836       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
837      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)') 
838      &  'Processor',myrank,': end reading molecular data.'
839 #endif
840       return
841       end
842 c--------------------------------------------------------------------------
843       logical function seq_comp(itypea,itypeb,length)
844       implicit none
845       integer length,itypea(length),itypeb(length)
846       integer i
847       do i=1,length
848         if (itypea(i).ne.itypeb(i)) then
849           seq_comp=.false.
850           return
851         endif
852       enddo
853       seq_comp=.true.
854       return
855       end
856 c-----------------------------------------------------------------------------
857       subroutine read_bridge
858 C Read information about disulfide bridges.
859       implicit real*8 (a-h,o-z)
860       include 'DIMENSIONS'
861 #ifdef MPI
862       include 'mpif.h'
863 #endif
864       include 'COMMON.IOUNITS'
865       include 'COMMON.GEO'
866       include 'COMMON.VAR'
867       include 'COMMON.INTERACT'
868       include 'COMMON.LOCAL'
869       include 'COMMON.NAMES'
870       include 'COMMON.CHAIN'
871       include 'COMMON.FFIELD'
872       include 'COMMON.SBRIDGE'
873       include 'COMMON.HEADER'
874       include 'COMMON.CONTROL'
875 c      include 'COMMON.DBASE'
876 c      include 'COMMON.THREAD'
877       include 'COMMON.TIME1'
878       include 'COMMON.SETUP'
879 C Read bridging residues.
880       read (inp,*) ns,(iss(i),i=1,ns)
881       print *,'ns=',ns
882       if(me.eq.king.or..not.out1file)
883      &  write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
884 C Check whether the specified bridging residues are cystines.
885       do i=1,ns
886         if (itype(iss(i)).ne.1) then
887           if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)') 
888      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
889      &   ' can form a disulfide bridge?!!!'
890           write (*,'(2a,i3,a)') 
891      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
892      &   ' can form a disulfide bridge?!!!'
893 #ifdef MPI
894          call MPI_Finalize(MPI_COMM_WORLD,ierror)
895          stop
896 #endif
897         endif
898       enddo
899 C Read preformed bridges.
900       if (ns.gt.0) then
901       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
902       write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
903       if (nss.gt.0) then
904         nhpb=nss
905 C Check if the residues involved in bridges are in the specified list of
906 C bridging residues.
907         do i=1,nss
908           do j=1,i-1
909             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
910      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
911               write (iout,'(a,i3,a)') 'Disulfide pair',i,
912      &      ' contains residues present in other pairs.'
913               write (*,'(a,i3,a)') 'Disulfide pair',i,
914      &      ' contains residues present in other pairs.'
915 #ifdef MPI
916               call MPI_Finalize(MPI_COMM_WORLD,ierror)
917               stop 
918 #endif
919             endif
920           enddo
921           do j=1,ns
922             if (ihpb(i).eq.iss(j)) goto 10
923           enddo
924           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
925    10     continue
926           do j=1,ns
927             if (jhpb(i).eq.iss(j)) goto 20
928           enddo
929           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
930    20     continue
931           dhpb(i)=dbr
932           forcon(i)=fbr
933         enddo
934         do i=1,nss
935           ihpb(i)=ihpb(i)+nres
936           jhpb(i)=jhpb(i)+nres
937         enddo
938       endif
939       endif
940       return
941       end
942 c----------------------------------------------------------------------------
943       subroutine read_x(kanal,*)
944       implicit real*8 (a-h,o-z)
945       include 'DIMENSIONS'
946       include 'COMMON.GEO'
947       include 'COMMON.VAR'
948       include 'COMMON.CHAIN'
949       include 'COMMON.IOUNITS'
950       include 'COMMON.CONTROL'
951       include 'COMMON.LOCAL'
952       include 'COMMON.INTERACT'
953 c Read coordinates from input
954 c
955       read(kanal,'(8f10.5)',end=10,err=10)
956      &  ((c(l,k),l=1,3),k=1,nres),
957      &  ((c(l,k+nres),l=1,3),k=nnt,nct)
958       do j=1,3
959         c(j,nres+1)=c(j,1)
960         c(j,2*nres)=c(j,nres)
961       enddo
962       call int_from_cart1(.false.)
963       do i=1,nres-1
964         do j=1,3
965           dc(j,i)=c(j,i+1)-c(j,i)
966           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
967         enddo
968       enddo
969       do i=nnt,nct
970         if (itype(i).ne.10) then
971           do j=1,3
972             dc(j,i+nres)=c(j,i+nres)-c(j,i)
973             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
974           enddo
975         endif
976       enddo
977
978       return
979    10 return1
980       end
981 c------------------------------------------------------------------------------
982       subroutine setup_var
983       implicit real*8 (a-h,o-z)
984       include 'DIMENSIONS'
985       include 'COMMON.IOUNITS'
986       include 'COMMON.GEO'
987       include 'COMMON.VAR'
988       include 'COMMON.INTERACT'
989       include 'COMMON.LOCAL'
990       include 'COMMON.NAMES'
991       include 'COMMON.CHAIN'
992       include 'COMMON.FFIELD'
993       include 'COMMON.SBRIDGE'
994       include 'COMMON.HEADER'
995       include 'COMMON.CONTROL'
996 c      include 'COMMON.DBASE'
997 c      include 'COMMON.THREAD'
998       include 'COMMON.TIME1'
999 C Set up variable list.
1000       ntheta=nres-2
1001       nphi=nres-3
1002       nvar=ntheta+nphi
1003       nside=0
1004       do i=2,nres-1
1005         if (itype(i).ne.10) then
1006           nside=nside+1
1007           ialph(i,1)=nvar+nside
1008           ialph(nside,2)=i
1009         endif
1010       enddo
1011       if (indphi.gt.0) then
1012         nvar=nphi
1013       else if (indback.gt.0) then
1014         nvar=nphi+ntheta
1015       else
1016         nvar=nvar+2*nside
1017       endif
1018 cd    write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1019       return
1020       end
1021 c----------------------------------------------------------------------------
1022       subroutine gen_dist_constr
1023 C Generate CA distance constraints.
1024       implicit real*8 (a-h,o-z)
1025       include 'DIMENSIONS'
1026       include 'COMMON.IOUNITS'
1027       include 'COMMON.GEO'
1028       include 'COMMON.VAR'
1029       include 'COMMON.INTERACT'
1030       include 'COMMON.LOCAL'
1031       include 'COMMON.NAMES'
1032       include 'COMMON.CHAIN'
1033       include 'COMMON.FFIELD'
1034       include 'COMMON.SBRIDGE'
1035       include 'COMMON.HEADER'
1036       include 'COMMON.CONTROL'
1037 c      include 'COMMON.DBASE'
1038 c      include 'COMMON.THREAD'
1039       include 'COMMON.TIME1'
1040       dimension itype_pdb(maxres)
1041       common /pizda/ itype_pdb
1042       character*2 iden
1043 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1044 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1045 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1046 cd     & ' nsup',nsup
1047       do i=nstart_sup,nstart_sup+nsup-1
1048 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1049 cd     &    ' seq_pdb', restyp(itype_pdb(i))
1050         do j=i+2,nstart_sup+nsup-1
1051           nhpb=nhpb+1
1052           ihpb(nhpb)=i+nstart_seq-nstart_sup
1053           jhpb(nhpb)=j+nstart_seq-nstart_sup
1054           forcon(nhpb)=weidis
1055           dhpb(nhpb)=dist(i,j)
1056         enddo
1057       enddo 
1058 cd      write (iout,'(a)') 'Distance constraints:' 
1059 cd      do i=nss+1,nhpb
1060 cd        ii=ihpb(i)
1061 cd        jj=jhpb(i)
1062 cd        iden='CA'
1063 cd        if (ii.gt.nres) then
1064 cd          iden='SC'
1065 cd          ii=ii-nres
1066 cd          jj=jj-nres
1067 cd        endif
1068 cd        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)') 
1069 cd     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1070 cd     &  dhpb(i),forcon(i)
1071 cd      enddo
1072       return
1073       end
1074 c----------------------------------------------------------------------------
1075       subroutine csaread
1076       implicit real*8 (a-h,o-z)
1077       include 'DIMENSIONS'
1078       include 'COMMON.IOUNITS'
1079       include 'COMMON.GEO'
1080       include 'COMMON.CSA'
1081       include 'COMMON.BANK'
1082       include 'COMMON.CONTROL'
1083       character*80 ucase
1084       character*620 mcmcard
1085       call card_concat(mcmcard)
1086
1087       call readi(mcmcard,'NCONF',nconf,50)
1088       call readi(mcmcard,'NADD',nadd,0)
1089       call readi(mcmcard,'JSTART',jstart,1)
1090       call readi(mcmcard,'JEND',jend,1)
1091       call readi(mcmcard,'NSTMAX',nstmax,500000)
1092       call readi(mcmcard,'N0',n0,1)
1093       call readi(mcmcard,'N1',n1,6)
1094       call readi(mcmcard,'N2',n2,4)
1095       call readi(mcmcard,'N3',n3,0)
1096       call readi(mcmcard,'N4',n4,0)
1097       call readi(mcmcard,'N5',n5,0)
1098       call readi(mcmcard,'N6',n6,10)
1099       call readi(mcmcard,'N7',n7,0)
1100       call readi(mcmcard,'N8',n8,0)
1101       call readi(mcmcard,'N9',n9,0)
1102       call readi(mcmcard,'N14',n14,0)
1103       call readi(mcmcard,'N15',n15,0)
1104       call readi(mcmcard,'N16',n16,0)
1105       call readi(mcmcard,'N17',n17,0)
1106       call readi(mcmcard,'N18',n18,0)
1107
1108       vdisulf=(index(mcmcard,'DYNSS').gt.0)
1109
1110       call readi(mcmcard,'NDIFF',ndiff,2)
1111       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1112       call readi(mcmcard,'IS1',is1,1)
1113       call readi(mcmcard,'IS2',is2,8)
1114       call readi(mcmcard,'NRAN0',nran0,4)
1115       call readi(mcmcard,'NRAN1',nran1,2)
1116       call readi(mcmcard,'IRR',irr,1)
1117       call readi(mcmcard,'NSEED',nseed,20)
1118       call readi(mcmcard,'NTOTAL',ntotal,10000)
1119       call reada(mcmcard,'CUT1',cut1,2.0d0)
1120       call reada(mcmcard,'CUT2',cut2,5.0d0)
1121       call reada(mcmcard,'ESTOP',estop,-300000.0d0)
1122       call readi(mcmcard,'ICMAX',icmax,1)
1123       call readi(mcmcard,'NBANKM',nbankm,400)
1124       call readi(mcmcard,'IUCUT',iucut,2)
1125       call readi(mcmcard,'IRESTART',irestart,0)
1126 c!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
1127       ntbankm=0
1128 c!bankt
1129       call reada(mcmcard,'DELE',dele,20.0d0)
1130       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1131       call readi(mcmcard,'IREF',iref,0)
1132       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1133       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1134       call readi(mcmcard,'NCONF_IN',nconf_in,0)
1135       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1136       write (iout,*) "NCONF_IN",nconf_in
1137       tm_score=(index(mcmcard,'TMSCORE').gt.0)
1138       if (tm_score) write (iout,*) "Using TM_Score instead of DIFF",
1139      &      " for torsional angles" 
1140       return
1141       end
1142
1143       subroutine read_minim
1144       implicit real*8 (a-h,o-z)
1145       include 'DIMENSIONS'
1146       include 'COMMON.MINIM'
1147       include 'COMMON.IOUNITS'
1148       character*80 ucase
1149       character*320 minimcard
1150       call card_concat(minimcard)
1151       call readi(minimcard,'MAXMIN',maxmin,2000)
1152       call readi(minimcard,'MAXFUN',maxfun,5000)
1153       call readi(minimcard,'MINMIN',minmin,maxmin)
1154       call readi(minimcard,'MINFUN',minfun,maxmin)
1155       call reada(minimcard,'TOLF',tolf,1.0D-2)
1156       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1157       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1158       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1159       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1160       write (iout,'(/80(1h*)/20x,a/80(1h*))') 
1161      &         'Options in energy minimization:'
1162       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1163      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1164      & 'MinMin:',MinMin,' MinFun:',MinFun,
1165      & ' TolF:',TolF,' RTolF:',RTolF
1166       return
1167       end
1168 c----------------------------------------------------------------------------
1169       subroutine read_angles(kanal,*)
1170       implicit real*8 (a-h,o-z)
1171       include 'DIMENSIONS'
1172       include 'COMMON.GEO'
1173       include 'COMMON.VAR'
1174       include 'COMMON.CHAIN'
1175       include 'COMMON.IOUNITS'
1176       include 'COMMON.CONTROL'
1177 c Read angles from input 
1178 c
1179        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1180        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1181        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1182        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1183
1184        do i=1,nres
1185 c 9/7/01 avoid 180 deg valence angle
1186         if (theta(i).gt.179.99d0) theta(i)=179.99d0
1187 c
1188         theta(i)=deg2rad*theta(i)
1189         phi(i)=deg2rad*phi(i)
1190         alph(i)=deg2rad*alph(i)
1191         omeg(i)=deg2rad*omeg(i)
1192        enddo
1193       return
1194    10 return1
1195       end
1196 c----------------------------------------------------------------------------
1197       subroutine reada(rekord,lancuch,wartosc,default)
1198       implicit none
1199       character*(*) rekord,lancuch
1200       double precision wartosc,default
1201       integer ilen,iread
1202       external ilen
1203       iread=index(rekord,lancuch)
1204       if (iread.eq.0) then
1205         wartosc=default 
1206         return
1207       endif   
1208       iread=iread+ilen(lancuch)+1
1209       read (rekord(iread:),*,err=10,end=10) wartosc
1210       return
1211   10  wartosc=default
1212       return
1213       end
1214 c----------------------------------------------------------------------------
1215       subroutine readi(rekord,lancuch,wartosc,default)
1216       implicit none
1217       character*(*) rekord,lancuch
1218       integer wartosc,default
1219       integer ilen,iread
1220       external ilen
1221       iread=index(rekord,lancuch)
1222       if (iread.eq.0) then
1223         wartosc=default 
1224         return
1225       endif   
1226       iread=iread+ilen(lancuch)+1
1227       read (rekord(iread:),*,err=10,end=10) wartosc
1228       return
1229   10  wartosc=default
1230       return
1231       end
1232 c----------------------------------------------------------------------------
1233       subroutine multreadi(rekord,lancuch,tablica,dim,default)
1234       implicit none
1235       integer dim,i
1236       integer tablica(dim),default
1237       character*(*) rekord,lancuch
1238       character*80 aux
1239       integer ilen,iread
1240       external ilen
1241       do i=1,dim
1242         tablica(i)=default
1243       enddo
1244       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1245       if (iread.eq.0) return
1246       iread=iread+ilen(lancuch)+1
1247       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1248    10 return
1249       end
1250 c----------------------------------------------------------------------------
1251       subroutine multreada(rekord,lancuch,tablica,dim,default)
1252       implicit none
1253       integer dim,i
1254       double precision tablica(dim),default
1255       character*(*) rekord,lancuch
1256       character*80 aux
1257       integer ilen,iread
1258       external ilen
1259       do i=1,dim
1260         tablica(i)=default
1261       enddo
1262       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1263       if (iread.eq.0) return
1264       iread=iread+ilen(lancuch)+1
1265       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1266    10 return
1267       end
1268 c----------------------------------------------------------------------------
1269       subroutine openunits
1270       implicit real*8 (a-h,o-z)
1271       include 'DIMENSIONS'    
1272 #ifdef MPI
1273       include 'mpif.h'
1274       character*16 form,nodename
1275       integer nodelen
1276 #endif
1277       include 'COMMON.SETUP'
1278       include 'COMMON.IOUNITS'
1279 c      include 'COMMON.MD'
1280       include 'COMMON.CONTROL'
1281       integer lenpre,lenpot,ilen,lentmp
1282       external ilen
1283       character*3 out1file_text,ucase
1284       character*3 ll
1285       external ucase
1286 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1287       call getenv_loc("PREFIX",prefix)
1288       pref_orig = prefix
1289       call getenv_loc("POT",pot)
1290       call getenv_loc("DIRTMP",tmpdir)
1291       call getenv_loc("CURDIR",curdir)
1292       call getenv_loc("OUT1FILE",out1file_text)
1293 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1294       out1file_text=ucase(out1file_text)
1295       if (out1file_text(1:1).eq."Y") then
1296         out1file=.true.
1297       else 
1298         out1file=fg_rank.gt.0
1299       endif
1300       lenpre=ilen(prefix)
1301       lenpot=ilen(pot)
1302       lentmp=ilen(tmpdir)
1303       if (lentmp.gt.0) then
1304           write (*,'(80(1h!))')
1305           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
1306           write (*,'(80(1h!))')
1307           write (*,*)"All output files will be on node /tmp directory." 
1308 #ifdef MPI
1309         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1310         if (me.eq.king) then
1311           write (*,*) "The master node is ",nodename
1312         else if (fg_rank.eq.0) then
1313           write (*,*) "I am the CG slave node ",nodename
1314         else 
1315           write (*,*) "I am the FG slave node ",nodename
1316         endif
1317 #endif
1318         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1319         lenpre = lentmp+lenpre+1
1320       endif
1321       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1322 C Get the names and open the input files
1323 #if defined(WINIFL) || defined(WINPGI)
1324       open(1,file=pref_orig(:ilen(pref_orig))//
1325      &  '.inp',status='old',readonly,shared)
1326        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1327 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1328 C Get parameter filenames and open the parameter files.
1329       call getenv_loc('BONDPAR',bondname)
1330       open (ibond,file=bondname,status='old',readonly,shared)
1331       call getenv_loc('THETPAR',thetname)
1332       open (ithep,file=thetname,status='old',readonly,shared)
1333 #ifndef CRYST_THETA
1334       call getenv_loc('THETPARPDB',thetname_pdb)
1335       open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
1336 #endif
1337       call getenv_loc('ROTPAR',rotname)
1338       open (irotam,file=rotname,status='old',readonly,shared)
1339 #ifndef CRYST_SC
1340       call getenv_loc('ROTPARPDB',rotname_pdb)
1341       open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
1342 #endif
1343       call getenv_loc('TORPAR',torname)
1344       open (itorp,file=torname,status='old',readonly,shared)
1345       call getenv_loc('TORDPAR',tordname)
1346       open (itordp,file=tordname,status='old',readonly,shared)
1347       call getenv_loc('FOURIER',fouriername)
1348       open (ifourier,file=fouriername,status='old',readonly,shared)
1349       call getenv_loc('ELEPAR',elename)
1350       open (ielep,file=elename,status='old',readonly,shared)
1351       call getenv_loc('SIDEPAR',sidename)
1352       open (isidep,file=sidename,status='old',readonly,shared)
1353 #elif (defined CRAY) || (defined AIX)
1354       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1355      &  action='read')
1356 c      print *,"Processor",myrank," opened file 1" 
1357       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1358 c      print *,"Processor",myrank," opened file 9" 
1359 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1360 C Get parameter filenames and open the parameter files.
1361       call getenv_loc('BONDPAR',bondname)
1362       open (ibond,file=bondname,status='old',action='read')
1363 c      print *,"Processor",myrank," opened file IBOND" 
1364       call getenv_loc('THETPAR',thetname)
1365       open (ithep,file=thetname,status='old',action='read')
1366 c      print *,"Processor",myrank," opened file ITHEP" 
1367 #ifndef CRYST_THETA
1368       call getenv_loc('THETPARPDB',thetname_pdb)
1369       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1370 #endif
1371       call getenv_loc('ROTPAR',rotname)
1372       open (irotam,file=rotname,status='old',action='read')
1373 c      print *,"Processor",myrank," opened file IROTAM" 
1374 #ifndef CRYST_SC
1375       call getenv_loc('ROTPARPDB',rotname_pdb)
1376       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
1377 #endif
1378       call getenv_loc('TORPAR',torname)
1379       open (itorp,file=torname,status='old',action='read')
1380 c      print *,"Processor",myrank," opened file ITORP" 
1381       call getenv_loc('TORDPAR',tordname)
1382       open (itordp,file=tordname,status='old',action='read')
1383 c      print *,"Processor",myrank," opened file ITORDP" 
1384       call getenv_loc('SCCORPAR',sccorname)
1385       open (isccor,file=sccorname,status='old',action='read')
1386 c      print *,"Processor",myrank," opened file ISCCOR" 
1387       call getenv_loc('FOURIER',fouriername)
1388       open (ifourier,file=fouriername,status='old',action='read')
1389 c      print *,"Processor",myrank," opened file IFOURIER" 
1390       call getenv_loc('ELEPAR',elename)
1391       open (ielep,file=elename,status='old',action='read')
1392 c      print *,"Processor",myrank," opened file IELEP" 
1393       call getenv_loc('SIDEPAR',sidename)
1394       open (isidep,file=sidename,status='old',action='read')
1395 c      print *,"Processor",myrank," opened file ISIDEP" 
1396 c      print *,"Processor",myrank," opened parameter files" 
1397 #elif (defined G77)
1398       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1399       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1400 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1401 C Get parameter filenames and open the parameter files.
1402       call getenv_loc('BONDPAR',bondname)
1403       open (ibond,file=bondname,status='old')
1404       call getenv_loc('THETPAR',thetname)
1405       open (ithep,file=thetname,status='old')
1406 #ifndef CRYST_THETA
1407       call getenv_loc('THETPARPDB',thetname_pdb)
1408       open (ithep_pdb,file=thetname_pdb,status='old')
1409 #endif
1410       call getenv_loc('ROTPAR',rotname)
1411       open (irotam,file=rotname,status='old')
1412 #ifndef CRYST_SC
1413       call getenv_loc('ROTPARPDB',rotname_pdb)
1414       open (irotam_pdb,file=rotname_pdb,status='old')
1415 #endif
1416       call getenv_loc('TORPAR',torname)
1417       open (itorp,file=torname,status='old')
1418       call getenv_loc('TORDPAR',tordname)
1419       open (itordp,file=tordname,status='old')
1420       call getenv_loc('SCCORPAR',sccorname)
1421       open (isccor,file=sccorname,status='old')
1422       call getenv_loc('FOURIER',fouriername)
1423       open (ifourier,file=fouriername,status='old')
1424       call getenv_loc('ELEPAR',elename)
1425       open (ielep,file=elename,status='old')
1426       call getenv_loc('SIDEPAR',sidename)
1427       open (isidep,file=sidename,status='old')
1428 #else
1429       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1430      &  readonly)
1431        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1432 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1433 C Get parameter filenames and open the parameter files.
1434       call getenv_loc('BONDPAR',bondname)
1435       open (ibond,file=bondname,status='old',readonly)
1436       call getenv_loc('THETPAR',thetname)
1437       open (ithep,file=thetname,status='old',readonly)
1438 #ifndef CRYST_THETA
1439       call getenv_loc('THETPARPDB',thetname_pdb)
1440       print *,"thetname_pdb ",thetname_pdb
1441       open (ithep_pdb,file=thetname_pdb,status='old',readonly)
1442       print *,ithep_pdb," opened"
1443 #endif
1444       call getenv_loc('ROTPAR',rotname)
1445       open (irotam,file=rotname,status='old',readonly)
1446 #ifndef CRYST_SC
1447       call getenv_loc('ROTPARPDB',rotname_pdb)
1448       open (irotam_pdb,file=rotname_pdb,status='old',readonly)
1449 #endif
1450       call getenv_loc('TORPAR',torname)
1451       open (itorp,file=torname,status='old',readonly)
1452       call getenv_loc('TORDPAR',tordname)
1453       open (itordp,file=tordname,status='old',readonly)
1454       call getenv_loc('SCCORPAR',sccorname)
1455       open (isccor,file=sccorname,status='old',readonly)
1456       call getenv_loc('FOURIER',fouriername)
1457       open (ifourier,file=fouriername,status='old',readonly)
1458       call getenv_loc('ELEPAR',elename)
1459       open (ielep,file=elename,status='old',readonly)
1460       call getenv_loc('SIDEPAR',sidename)
1461       open (isidep,file=sidename,status='old',readonly)
1462 #endif
1463 #ifndef OLDSCP
1464 C
1465 C 8/9/01 In the newest version SCp interaction constants are read from a file
1466 C Use -DOLDSCP to use hard-coded constants instead.
1467 C
1468       call getenv_loc('SCPPAR',scpname)
1469 #if defined(WINIFL) || defined(WINPGI)
1470       open (iscpp,file=scpname,status='old',readonly,shared)
1471 #elif (defined CRAY)  || (defined AIX)
1472       open (iscpp,file=scpname,status='old',action='read')
1473 #elif (defined G77)
1474       open (iscpp,file=scpname,status='old')
1475 #else
1476       open (iscpp,file=scpname,status='old',readonly)
1477 #endif
1478 #endif
1479       call getenv_loc('PATTERN',patname)
1480 #if defined(WINIFL) || defined(WINPGI)
1481       open (icbase,file=patname,status='old',readonly,shared)
1482 #elif (defined CRAY)  || (defined AIX)
1483       open (icbase,file=patname,status='old',action='read')
1484 #elif (defined G77)
1485       open (icbase,file=patname,status='old')
1486 #else
1487       open (icbase,file=patname,status='old',readonly)
1488 #endif
1489 #ifdef MPI
1490 C Open output file only for CG processes
1491 c      print *,"Processor",myrank," fg_rank",fg_rank
1492       if (fg_rank.eq.0) then
1493
1494       if (nodes.eq.1) then
1495         npos=3
1496       else
1497         npos = dlog10(dfloat(nodes-1))+1
1498       endif
1499       if (npos.lt.3) npos=3
1500       write (liczba,'(i1)') npos
1501       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
1502      &  //')'
1503       write (liczba,form) me
1504       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
1505      &  liczba(:ilen(liczba))
1506       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1507      &  //'.int'
1508       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1509      &  //'.pdb'
1510       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
1511      &  liczba(:ilen(liczba))//'.mol2'
1512       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
1513      &  liczba(:ilen(liczba))//'.stat'
1514       if (lentmp.gt.0)
1515      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
1516      &      //liczba(:ilen(liczba))//'.stat')
1517       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
1518      &  //'.rst'
1519 c      if(usampl) then
1520 c          qname=prefix(:lenpre)//'_'//pot(:lenpot)//
1521 c     & liczba(:ilen(liczba))//'.const'
1522 c      endif 
1523
1524       endif
1525 #else
1526       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
1527       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
1528       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
1529       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
1530       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
1531       if (lentmp.gt.0)
1532      & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)//
1533      &    '.stat')
1534       rest2name=prefix(:ilen(prefix))//'.rst'
1535 c      if(usampl) then 
1536 c         qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
1537 c      endif 
1538 #endif
1539 #if defined(AIX) || defined(PGI)
1540       if (me.eq.king .or. .not. out1file) 
1541      &   open(iout,file=outname,status='unknown')
1542 #ifdef DEBUG
1543       if (fg_rank.gt.0) then
1544         write (liczba,'(i3.3)') myrank/nfgtasks
1545         write (ll,'(bz,i3.3)') fg_rank
1546         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1547      &   status='unknown')
1548       endif
1549 #endif
1550       if(me.eq.king) then
1551        open(igeom,file=intname,status='unknown',position='append')
1552        open(ipdb,file=pdbname,status='unknown')
1553        open(imol2,file=mol2name,status='unknown')
1554        open(istat,file=statname,status='unknown',position='append')
1555       else
1556 c1out       open(iout,file=outname,status='unknown')
1557       endif
1558 #else
1559       if (me.eq.king .or. .not.out1file)
1560      &    open(iout,file=outname,status='unknown')
1561 #ifdef DEBUG
1562       if (fg_rank.gt.0) then
1563         write (liczba,'(i3.3)') myrank/nfgtasks
1564         write (ll,'(bz,i3.3)') fg_rank
1565         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1566      &   status='unknown')
1567       endif
1568 #endif
1569       if(me.eq.king) then
1570        open(igeom,file=intname,status='unknown',access='append')
1571        open(ipdb,file=pdbname,status='unknown')
1572        open(imol2,file=mol2name,status='unknown')
1573        open(istat,file=statname,status='unknown',access='append')
1574       else
1575 c1out       open(iout,file=outname,status='unknown')
1576       endif
1577 #endif
1578       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
1579       csa_seed=prefix(:lenpre)//'.CSA.seed'
1580       csa_history=prefix(:lenpre)//'.CSA.history'
1581       csa_bank=prefix(:lenpre)//'.CSA.bank'
1582       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
1583       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
1584       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
1585 c!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
1586       csa_int=prefix(:lenpre)//'.int'
1587       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
1588       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
1589       csa_in=prefix(:lenpre)//'.CSA.in'
1590 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
1591 C Write file names
1592       if (me.eq.king)then
1593       write (iout,'(80(1h-))')
1594       write (iout,'(30x,a)') "FILE ASSIGNMENT"
1595       write (iout,'(80(1h-))')
1596       write (iout,*) "Input file                      : ",
1597      &  pref_orig(:ilen(pref_orig))//'.inp'
1598       write (iout,*) "Output file                     : ",
1599      &  outname(:ilen(outname))
1600       write (iout,*)
1601       write (iout,*) "Sidechain potential file        : ",
1602      &  sidename(:ilen(sidename))
1603 #ifndef OLDSCP
1604       write (iout,*) "SCp potential file              : ",
1605      &  scpname(:ilen(scpname))
1606 #endif
1607       write (iout,*) "Electrostatic potential file    : ",
1608      &  elename(:ilen(elename))
1609       write (iout,*) "Cumulant coefficient file       : ",
1610      &  fouriername(:ilen(fouriername))
1611       write (iout,*) "Torsional parameter file        : ",
1612      &  torname(:ilen(torname))
1613       write (iout,*) "Double torsional parameter file : ",
1614      &  tordname(:ilen(tordname))
1615       write (iout,*) "SCCOR parameter file : ",
1616      &  sccorname(:ilen(sccorname))
1617       write (iout,*) "Bond & inertia constant file    : ",
1618      &  bondname(:ilen(bondname))
1619       write (iout,*) "Bending parameter file          : ",
1620      &  thetname(:ilen(thetname))
1621       write (iout,*) "Rotamer parameter file          : ",
1622      &  rotname(:ilen(rotname))
1623       write (iout,*) "Threading database              : ",
1624      &  patname(:ilen(patname))
1625       if (lentmp.ne.0) 
1626      &write (iout,*)" DIRTMP                          : ",
1627      &  tmpdir(:lentmp)
1628       write (iout,'(80(1h-))')
1629       endif
1630       return
1631       end
1632 c----------------------------------------------------------------------------
1633       subroutine card_concat(card)
1634       implicit real*8 (a-h,o-z)
1635       include 'DIMENSIONS'
1636       include 'COMMON.IOUNITS'
1637       character*(*) card
1638       character*80 karta,ucase
1639       external ilen
1640       read (inp,'(a)') karta
1641       karta=ucase(karta)
1642       card=' '
1643       do while (karta(80:80).eq.'&')
1644         card=card(:ilen(card)+1)//karta(:79)
1645         read (inp,'(a)') karta
1646         karta=ucase(karta)
1647       enddo
1648       card=card(:ilen(card)+1)//karta
1649       return
1650       end
1651
1652       subroutine read_dist_constr
1653       implicit real*8 (a-h,o-z)
1654       include 'DIMENSIONS'
1655 #ifdef MPI
1656       include 'mpif.h'
1657 #endif
1658       include 'COMMON.SETUP'
1659       include 'COMMON.CONTROL'
1660       include 'COMMON.CHAIN'
1661       include 'COMMON.IOUNITS'
1662       include 'COMMON.SBRIDGE'
1663       integer ifrag_(2,100),ipair_(2,100)
1664       double precision wfrag_(100),wpair_(100)
1665       character*500 controlcard
1666       write (iout,*) "Calling read_dist_constr"
1667       write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
1668       call flush(iout)
1669       call card_concat(controlcard)
1670       call readi(controlcard,"NFRAG",nfrag_,0)
1671       call readi(controlcard,"NPAIR",npair_,0)
1672       call readi(controlcard,"NDIST",ndist_,0)
1673       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
1674       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
1675       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
1676       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
1677       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
1678 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
1679 c      write (iout,*) "IFRAG"
1680 c      do i=1,nfrag_
1681 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1682 c      enddo
1683 c      write (iout,*) "IPAIR"
1684 c      do i=1,npair_
1685 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
1686 c      enddo
1687       call flush(iout)
1688       do i=1,nfrag_
1689         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
1690         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
1691      &    ifrag_(2,i)=nstart_sup+nsup-1
1692 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1693         call flush(iout)
1694         if (wfrag_(i).gt.0.0d0) then
1695         do j=ifrag_(1,i),ifrag_(2,i)-1
1696           do k=j+1,ifrag_(2,i)
1697             write (iout,*) "j",j," k",k
1698             ddjk=dist(j,k)
1699             if (constr_dist.eq.1) then
1700             nhpb=nhpb+1
1701             ihpb(nhpb)=j
1702             jhpb(nhpb)=k
1703               dhpb(nhpb)=ddjk
1704             forcon(nhpb)=wfrag_(i) 
1705             else if (constr_dist.eq.2) then
1706               if (ddjk.le.dist_cut) then
1707                 nhpb=nhpb+1
1708                 ihpb(nhpb)=j
1709                 jhpb(nhpb)=k
1710                 dhpb(nhpb)=ddjk
1711                 forcon(nhpb)=wfrag_(i) 
1712               endif
1713             else
1714               nhpb=nhpb+1
1715               ihpb(nhpb)=j
1716               jhpb(nhpb)=k
1717               dhpb(nhpb)=ddjk
1718               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1719             endif
1720 #ifdef MPI
1721             if (.not.out1file .or. me.eq.king) 
1722      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1723      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1724 #else
1725             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1726      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1727 #endif
1728           enddo
1729         enddo
1730         endif
1731       enddo
1732       do i=1,npair_
1733         if (wpair_(i).gt.0.0d0) then
1734         ii = ipair_(1,i)
1735         jj = ipair_(2,i)
1736         if (ii.gt.jj) then
1737           itemp=ii
1738           ii=jj
1739           jj=itemp
1740         endif
1741         do j=ifrag_(1,ii),ifrag_(2,ii)
1742           do k=ifrag_(1,jj),ifrag_(2,jj)
1743             nhpb=nhpb+1
1744             ihpb(nhpb)=j
1745             jhpb(nhpb)=k
1746             forcon(nhpb)=wpair_(i)
1747             dhpb(nhpb)=dist(j,k)
1748 #ifdef MPI
1749             if (.not.out1file .or. me.eq.king)
1750      &      write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1751      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1752 #else
1753             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1754      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1755 #endif
1756           enddo
1757         enddo
1758         endif
1759       enddo 
1760       do i=1,ndist_
1761         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
1762         if (forcon(nhpb+1).gt.0.0d0) then
1763           nhpb=nhpb+1
1764           dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1765 #ifdef MPI
1766           if (.not.out1file .or. me.eq.king)
1767      &    write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1768      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1769 #else
1770           write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1771      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1772 #endif
1773         endif
1774       enddo
1775       call flush(iout)
1776       return
1777       end
1778 c-------------------------------------------------------------------------------
1779 #ifdef WINIFL
1780       subroutine flush(iu)
1781       return
1782       end
1783 #endif
1784 #ifdef AIX
1785       subroutine flush(iu)
1786       call flush_(iu)
1787       return
1788       end
1789 #endif
1790 c------------------------------------------------------------------------------
1791       subroutine copy_to_tmp(source)
1792       include "DIMENSIONS"
1793       include "COMMON.IOUNITS"
1794       character*(*) source
1795       character* 256 tmpfile
1796       integer ilen
1797       external ilen
1798       logical ex
1799       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
1800       inquire(file=tmpfile,exist=ex)
1801       if (ex) then
1802         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
1803      &   " to temporary directory..."
1804         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
1805         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
1806       endif
1807       return
1808       end
1809 c------------------------------------------------------------------------------
1810       subroutine move_from_tmp(source)
1811       include "DIMENSIONS"
1812       include "COMMON.IOUNITS"
1813       character*(*) source
1814       integer ilen
1815       external ilen
1816       write (*,*) "Moving ",source(:ilen(source)),
1817      & " from temporary directory to working directory"
1818       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
1819       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
1820       return
1821       end
1822 c------------------------------------------------------------------------------
1823       subroutine random_init(seed)
1824 C
1825 C Initialize random number generator
1826 C
1827       implicit real*8 (a-h,o-z)
1828       include 'DIMENSIONS'
1829 #ifdef AMD64
1830       integer*8 iseedi8
1831 #endif
1832 #ifdef MPI
1833       include 'mpif.h'
1834       logical OKRandom, prng_restart
1835       real*8  r1
1836       integer iseed_array(4)
1837 #endif
1838       include 'COMMON.IOUNITS'
1839       include 'COMMON.TIME1'
1840 c      include 'COMMON.THREAD'
1841       include 'COMMON.SBRIDGE'
1842       include 'COMMON.CONTROL'
1843       include 'COMMON.MCM'
1844 c      include 'COMMON.MAP'
1845       include 'COMMON.HEADER'
1846 c      include 'COMMON.CSA'
1847       include 'COMMON.CHAIN'
1848 c      include 'COMMON.MUCA'
1849 c      include 'COMMON.MD'
1850       include 'COMMON.FFIELD'
1851       include 'COMMON.SETUP'
1852       iseed=-dint(dabs(seed))
1853       if (iseed.eq.0) then
1854         write (iout,'(/80(1h*)/20x,a/80(1h*))') 
1855      &    'Random seed undefined. The program will stop.'
1856         write (*,'(/80(1h*)/20x,a/80(1h*))') 
1857      &    'Random seed undefined. The program will stop.'
1858 #ifdef MPI
1859         call mpi_finalize(mpi_comm_world,ierr)
1860 #endif
1861         stop 'Bad random seed.'
1862       endif
1863 #ifdef MPI
1864       if (fg_rank.eq.0) then
1865       seed=seed*(me+1)+1
1866 #ifdef AMD64
1867       iseedi8=dint(seed)
1868       if(me.eq.king .or. .not. out1file)
1869      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1870       write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1871       OKRandom = prng_restart(me,iseedi8)
1872 #else
1873       do i=1,4
1874        tmp=65536.0d0**(4-i)
1875        iseed_array(i) = dint(seed/tmp)
1876        seed=seed-iseed_array(i)*tmp
1877       enddo
1878       if(me.eq.king .or. .not. out1file)
1879      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
1880      &                 (iseed_array(i),i=1,4)
1881       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
1882      &                 (iseed_array(i),i=1,4)
1883       OKRandom = prng_restart(me,iseed_array)
1884 #endif
1885       if (OKRandom) then
1886         r1=ran_number(0.0D0,1.0D0)
1887         if(me.eq.king .or. .not. out1file)
1888      &   write (iout,*) 'ran_num',r1
1889         if (r1.lt.0.0d0) OKRandom=.false.
1890       endif
1891       if (.not.OKRandom) then
1892         write (iout,*) 'PRNG IS NOT WORKING!!!'
1893         print *,'PRNG IS NOT WORKING!!!'
1894         if (me.eq.0) then 
1895          call flush(iout)
1896          call mpi_abort(mpi_comm_world,error_msg,ierr)
1897          stop
1898         else
1899          write (iout,*) 'too many processors for parallel prng'
1900          write (*,*) 'too many processors for parallel prng'
1901          call flush(iout)
1902          stop
1903         endif
1904       endif
1905       endif
1906 #else
1907       call vrndst(iseed)
1908       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
1909 #endif
1910       return
1911       end