energy_dec for esccor and some other debug printout
[unres.git] / source / unres / src_CSA / 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       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 C Following 2 lines for diagnostics; comment out if not needed
459          write (iout,*) "Before sideadd"
460          call intout
461          if(me.eq.king.or..not.out1file)
462      &    write(iout,*)'Adding sidechains'
463          maxsi=1000
464          do i=2,nres-1
465           iti=itype(i)
466           if (iti.ne.10) then
467             nsi=0
468             fail=.true.
469             do while (fail.and.nsi.le.maxsi)
470 c              call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
471               nsi=nsi+1
472             enddo
473             if(fail) write(iout,*)'Adding sidechain failed for res ',
474      &              i,' after ',nsi,' trials'
475           endif
476          enddo
477         endif  
478 C 10/03/12 Adam: Recalculate coordinates with new side chain positions
479         call chainbuild
480 C Following 2 lines for diagnostics; comment out if not needed
481         write (iout,*) "After sideadd"
482         call intout
483       endif
484
485       if (indpdb.eq.0) then
486 C Read sequence if not taken from the pdb file.
487         read (inp,*) nres
488 c        print *,'nres=',nres
489         if (iscode.gt.0) then
490           read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
491         else
492           read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
493         endif
494 C Convert sequence to numeric code
495         do i=1,nres
496           itype(i)=rescode(i,sequence(i),iscode)
497         enddo
498 C Assign initial virtual bond lengths
499         do i=2,nres
500           vbld(i)=vbl
501           vbld_inv(i)=vblinv
502         enddo
503         do i=2,nres-1
504           vbld(i+nres)=dsc(itype(i))
505           vbld_inv(i+nres)=dsc_inv(itype(i))
506 c          write (iout,*) "i",i," itype",itype(i),
507 c     &      " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
508         enddo
509       endif 
510 c      print *,nres
511 c      print '(20i4)',(itype(i),i=1,nres)
512       do i=1,nres
513 #ifdef PROCOR
514         if (itype(i).eq.21 .or. itype(i+1).eq.21) then
515 #else
516         if (itype(i).eq.21) then
517 #endif
518           itel(i)=0
519 #ifdef PROCOR
520         else if (itype(i+1).ne.20) then
521 #else
522         else if (itype(i).ne.20) then
523 #endif
524           itel(i)=1
525         else
526           itel(i)=2
527         endif  
528       enddo
529       if(me.eq.king.or..not.out1file)then
530        write (iout,*) "ITEL"
531        do i=1,nres-1
532          write (iout,*) i,itype(i),itel(i)
533        enddo
534        print *,'Call Read_Bridge.'
535       endif
536       call read_bridge
537 C 8/13/98 Set limits to generating the dihedral angles
538       do i=1,nres
539         phibound(1,i)=-pi
540         phibound(2,i)=pi
541       enddo
542       read (inp,*) ndih_constr
543       if (ndih_constr.gt.0) then
544         read (inp,*) ftors
545         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
546         if(me.eq.king.or..not.out1file)then
547          write (iout,*) 
548      &   'There are',ndih_constr,' constraints on phi angles.'
549          do i=1,ndih_constr
550           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
551          enddo
552         endif
553         do i=1,ndih_constr
554           phi0(i)=deg2rad*phi0(i)
555           drange(i)=deg2rad*drange(i)
556         enddo
557         if(me.eq.king.or..not.out1file)
558      &   write (iout,*) 'FTORS',ftors
559         do i=1,ndih_constr
560           ii = idih_constr(i)
561           phibound(1,ii) = phi0(i)-drange(i)
562           phibound(2,ii) = phi0(i)+drange(i)
563         enddo 
564       endif
565       nnt=1
566 #ifdef MPI
567       if (me.eq.king) then
568 #endif
569        write (iout,'(a)') 'Boundaries in phi angle sampling:'
570        do i=1,nres
571          write (iout,'(a3,i5,2f10.1)') 
572      &   restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
573        enddo
574 #ifdef MP
575       endif
576 #endif
577       nct=nres
578 cd      print *,'NNT=',NNT,' NCT=',NCT
579       if (itype(1).eq.21) nnt=2
580       if (itype(nres).eq.21) nct=nct-1
581
582 C     Juyong:READ init_vars
583 C     Initialize variables!
584 C     Juyong:READ read_info
585 C     READ fragment information!!
586 C     both routines should be in dfa.F file!!
587
588       if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
589      &            wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
590        call init_dfa_vars
591        print*, 'init_dfa_vars finished!'
592        call read_dfa_info
593        print*, 'read_dfa_info finished!'
594       endif
595 C
596 C
597
598
599       if (pdbref) then
600         if(me.eq.king.or..not.out1file)
601      &   write (iout,'(a,i3)') 'nsup=',nsup
602         nstart_seq=nnt
603         if (nsup.le.(nct-nnt+1)) then
604           do i=0,nct-nnt+1-nsup
605             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
606               nstart_seq=nnt+i
607               goto 111
608             endif
609           enddo
610           write (iout,'(a)') 
611      &            'Error - sequences to be superposed do not match.'
612           stop
613         else
614           do i=0,nsup-(nct-nnt+1)
615             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1)) 
616      &      then
617               nstart_sup=nstart_sup+i
618               nsup=nct-nnt+1
619               goto 111
620             endif
621           enddo 
622           write (iout,'(a)') 
623      &            'Error - sequences to be superposed do not match.'
624         endif
625   111   continue
626         if (nsup.eq.0) nsup=nct-nnt
627         if (nstart_sup.eq.0) nstart_sup=nnt
628         if (nstart_seq.eq.0) nstart_seq=nnt
629         if(me.eq.king.or..not.out1file)  
630      &   write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
631      &                 ' nstart_seq=',nstart_seq
632       endif
633 c--- Zscore rms -------
634       if (nz_start.eq.0) nz_start=nnt
635       if (nz_end.eq.0 .and. nsup.gt.0) then
636         nz_end=nnt+nsup-1
637       else if (nz_end.eq.0) then
638         nz_end=nct
639       endif
640       if(me.eq.king.or..not.out1file)then
641        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
642        write (iout,*) 'IZ_SC=',iz_sc
643       endif
644 c----------------------
645       call init_int_table
646       if (refstr) then
647         if (.not.pdbref) then
648           call read_angles(inp,*38)
649           goto 39
650    38     write (iout,'(a)') 'Error reading reference structure.'
651 #ifdef MPI
652           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
653           stop 'Error reading reference structure'
654 #endif
655    39     call chainbuild
656           call setup_var
657 czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
658           nstart_sup=nnt
659           nstart_seq=nnt
660           nsup=nct-nnt+1
661           do i=1,2*nres
662             do j=1,3
663               cref(j,i)=c(j,i)
664             enddo
665           enddo
666           call contact(.true.,ncont_ref,icont_ref,co)
667         endif
668         if(me.eq.king.or..not.out1file)
669      &   write (iout,*) 'Contact order:',co
670         if (pdbref) then
671         if(me.eq.king.or..not.out1file)
672      &   write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
673         do i=1,ncont_ref
674           do j=1,2
675             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
676           enddo
677           if(me.eq.king.or..not.out1file)
678      &     write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
679      &     icont_ref(1,i),' ',
680      &     restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
681         enddo
682         endif
683       endif
684 c        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
685       if (constr_dist.gt.0) then
686         call read_dist_constr
687       endif
688       if (nhpb.gt.0) call hpb_partition
689 c      write (iout,*) "After read_dist_constr nhpb",nhpb
690 c      call flush(iout)
691       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
692      &    .and. modecalc.ne.8 .and. modecalc.ne.9 .and. 
693      &    modecalc.ne.10) then
694 C If input structure hasn't been supplied from the PDB file read or generate
695 C initial geometry.
696         if (iranconf.eq.0 .and. .not. extconf) then
697           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
698      &     write (iout,'(a)') 'Initial geometry will be read in.'
699           if (read_cart) then
700             read(inp,'(8f10.5)',end=36,err=36)
701      &       ((c(l,k),l=1,3),k=1,nres),
702      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
703             call int_from_cart1(.false.)
704             do i=1,nres-1
705               do j=1,3
706                 dc(j,i)=c(j,i+1)-c(j,i)
707                 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
708               enddo
709             enddo
710             do i=nnt,nct
711               if (itype(i).ne.10) then
712                 do j=1,3
713                   dc(j,i+nres)=c(j,i+nres)-c(j,i) 
714                   dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
715                 enddo
716               endif
717             enddo
718             return
719           else
720             call read_angles(inp,*36)
721           endif
722           goto 37
723    36     write (iout,'(a)') 'Error reading angle file.'
724 #ifdef MPI
725           call mpi_finalize( MPI_COMM_WORLD,IERR )
726 #endif
727           stop 'Error reading angle file.'
728    37     continue 
729         else if (extconf) then
730          if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
731      &    write (iout,'(a)') 'Extended chain initial geometry.'
732          do i=3,nres
733           theta(i)=90d0*deg2rad
734          enddo
735          do i=4,nres
736           phi(i)=180d0*deg2rad
737          enddo
738          do i=2,nres-1
739           alph(i)=110d0*deg2rad
740          enddo
741          do i=2,nres-1
742           omeg(i)=-120d0*deg2rad
743          enddo
744         else
745           if(me.eq.king.or..not.out1file)
746      &     write (iout,'(a)') 'Random-generated initial geometry.'
747
748
749 #ifdef MPI
750           if (me.eq.king  .or. fg_rank.eq.0 .and. (
751      &           modecalc.eq.12 .or. modecalc.eq.14) ) then  
752 #endif
753             do itrial=1,100
754               itmp=1
755 c              call gen_rand_conf(itmp,*30)
756               goto 40
757    30         write (iout,*) 'Failed to generate random conformation',
758      &          ', itrial=',itrial
759               write (*,*) 'Processor:',me,
760      &          ' Failed to generate random conformation',
761      &          ' itrial=',itrial
762               call intout
763
764 #ifdef AIX
765               call flush_(iout)
766 #else
767               call flush(iout)
768 #endif
769             enddo
770             write (iout,'(a,i3,a)') 'Processor:',me,
771      &        ' error in generating random conformation.'
772             write (*,'(a,i3,a)') 'Processor:',me,
773      &        ' error in generating random conformation.'
774             call flush(iout)
775 #ifdef MPI
776             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)            
777    40       continue
778           endif
779 #else
780           do itrial=1,100
781             itmp=1
782 c            call gen_rand_conf(itmp,*31)
783             goto 40
784    31       write (iout,*) 'Failed to generate random conformation',
785      &        ', itrial=',itrial
786             write (*,*) 'Failed to generate random conformation',
787      &        ', itrial=',itrial
788           enddo
789           write (iout,'(a,i3,a)') 'Processor:',me,
790      &      ' error in generating random conformation.'
791           write (*,'(a,i3,a)') 'Processor:',me,
792      &      ' error in generating random conformation.'
793           stop
794    40     continue
795 #endif
796         endif
797       elseif (modecalc.eq.4) then
798         read (inp,'(a)') intinname
799         open (intin,file=intinname,status='old',err=333)
800         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
801      &  write (iout,'(a)') 'intinname',intinname
802         write (*,'(a)') 'Processor',myrank,' intinname',intinname
803         goto 334
804   333   write (iout,'(2a)') 'Error opening angle file ',intinname
805 #ifdef MPI 
806         call MPI_Finalize(MPI_COMM_WORLD,IERR)
807 #endif   
808         stop 'Error opening angle file.' 
809   334   continue
810
811       endif 
812 C Generate distance constraints, if the PDB structure is to be regularized. 
813       call setup_var
814       if (me.eq.king .or. .not. out1file)
815      & call intout
816       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
817         write (iout,'(/a,i3,a)') 
818      &  'The chain contains',ns,' disulfide-bridging cysteines.'
819         write (iout,'(20i4)') (iss(i),i=1,ns)
820         write (iout,'(/a/)') 'Pre-formed links are:' 
821         do i=1,nss
822           i1=ihpb(i)-nres
823           i2=jhpb(i)-nres
824           it1=itype(i1)
825           it2=itype(i2)
826           if (me.eq.king.or..not.out1file)
827      &    write (iout,'(2a,i3,3a,i3,a,3f10.3)')
828      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
829      &    ebr,forcon(i)
830         enddo
831         write (iout,'(a)')
832       endif
833 c       if (i2ndstr.gt.0) call secstrp2dihc
834 c      call geom_to_var(nvar,x)
835 c      call etotal(energia(0))
836 c      call enerprint(energia(0))
837 c      call briefout(0,etot)
838 c      stop
839 cd    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
840 cd    write (iout,'(a)') 'Variable list:'
841 cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
842 #ifdef MPI
843       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
844      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)') 
845      &  'Processor',myrank,': end reading molecular data.'
846 #endif
847       return
848       end
849 c--------------------------------------------------------------------------
850       logical function seq_comp(itypea,itypeb,length)
851       implicit none
852       integer length,itypea(length),itypeb(length)
853       integer i
854       do i=1,length
855         if (itypea(i).ne.itypeb(i)) then
856           seq_comp=.false.
857           return
858         endif
859       enddo
860       seq_comp=.true.
861       return
862       end
863 c-----------------------------------------------------------------------------
864       subroutine read_bridge
865 C Read information about disulfide bridges.
866       implicit real*8 (a-h,o-z)
867       include 'DIMENSIONS'
868 #ifdef MPI
869       include 'mpif.h'
870 #endif
871       include 'COMMON.IOUNITS'
872       include 'COMMON.GEO'
873       include 'COMMON.VAR'
874       include 'COMMON.INTERACT'
875       include 'COMMON.LOCAL'
876       include 'COMMON.NAMES'
877       include 'COMMON.CHAIN'
878       include 'COMMON.FFIELD'
879       include 'COMMON.SBRIDGE'
880       include 'COMMON.HEADER'
881       include 'COMMON.CONTROL'
882 c      include 'COMMON.DBASE'
883 c      include 'COMMON.THREAD'
884       include 'COMMON.TIME1'
885       include 'COMMON.SETUP'
886 C Read bridging residues.
887       read (inp,*) ns,(iss(i),i=1,ns)
888       print *,'ns=',ns
889       if(me.eq.king.or..not.out1file)
890      &  write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
891 C Check whether the specified bridging residues are cystines.
892       do i=1,ns
893         if (itype(iss(i)).ne.1) then
894           if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)') 
895      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
896      &   ' can form a disulfide bridge?!!!'
897           write (*,'(2a,i3,a)') 
898      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
899      &   ' can form a disulfide bridge?!!!'
900 #ifdef MPI
901          call MPI_Finalize(MPI_COMM_WORLD,ierror)
902          stop
903 #endif
904         endif
905       enddo
906 C Read preformed bridges.
907       if (ns.gt.0) then
908       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
909       write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
910       if (nss.gt.0) then
911         nhpb=nss
912 C Check if the residues involved in bridges are in the specified list of
913 C bridging residues.
914         do i=1,nss
915           do j=1,i-1
916             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
917      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
918               write (iout,'(a,i3,a)') 'Disulfide pair',i,
919      &      ' contains residues present in other pairs.'
920               write (*,'(a,i3,a)') 'Disulfide pair',i,
921      &      ' contains residues present in other pairs.'
922 #ifdef MPI
923               call MPI_Finalize(MPI_COMM_WORLD,ierror)
924               stop 
925 #endif
926             endif
927           enddo
928           do j=1,ns
929             if (ihpb(i).eq.iss(j)) goto 10
930           enddo
931           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
932    10     continue
933           do j=1,ns
934             if (jhpb(i).eq.iss(j)) goto 20
935           enddo
936           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
937    20     continue
938           dhpb(i)=dbr
939           forcon(i)=fbr
940         enddo
941         do i=1,nss
942           ihpb(i)=ihpb(i)+nres
943           jhpb(i)=jhpb(i)+nres
944         enddo
945       endif
946       endif
947       return
948       end
949 c----------------------------------------------------------------------------
950       subroutine read_x(kanal,*)
951       implicit real*8 (a-h,o-z)
952       include 'DIMENSIONS'
953       include 'COMMON.GEO'
954       include 'COMMON.VAR'
955       include 'COMMON.CHAIN'
956       include 'COMMON.IOUNITS'
957       include 'COMMON.CONTROL'
958       include 'COMMON.LOCAL'
959       include 'COMMON.INTERACT'
960 c Read coordinates from input
961 c
962       read(kanal,'(8f10.5)',end=10,err=10)
963      &  ((c(l,k),l=1,3),k=1,nres),
964      &  ((c(l,k+nres),l=1,3),k=nnt,nct)
965       do j=1,3
966         c(j,nres+1)=c(j,1)
967         c(j,2*nres)=c(j,nres)
968       enddo
969       call int_from_cart1(.false.)
970       do i=1,nres-1
971         do j=1,3
972           dc(j,i)=c(j,i+1)-c(j,i)
973           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
974         enddo
975       enddo
976       do i=nnt,nct
977         if (itype(i).ne.10) then
978           do j=1,3
979             dc(j,i+nres)=c(j,i+nres)-c(j,i)
980             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
981           enddo
982         endif
983       enddo
984
985       return
986    10 return1
987       end
988 c------------------------------------------------------------------------------
989       subroutine setup_var
990       implicit real*8 (a-h,o-z)
991       include 'DIMENSIONS'
992       include 'COMMON.IOUNITS'
993       include 'COMMON.GEO'
994       include 'COMMON.VAR'
995       include 'COMMON.INTERACT'
996       include 'COMMON.LOCAL'
997       include 'COMMON.NAMES'
998       include 'COMMON.CHAIN'
999       include 'COMMON.FFIELD'
1000       include 'COMMON.SBRIDGE'
1001       include 'COMMON.HEADER'
1002       include 'COMMON.CONTROL'
1003 c      include 'COMMON.DBASE'
1004 c      include 'COMMON.THREAD'
1005       include 'COMMON.TIME1'
1006 C Set up variable list.
1007       ntheta=nres-2
1008       nphi=nres-3
1009       nvar=ntheta+nphi
1010       nside=0
1011       do i=2,nres-1
1012         if (itype(i).ne.10) then
1013           nside=nside+1
1014           ialph(i,1)=nvar+nside
1015           ialph(nside,2)=i
1016         endif
1017       enddo
1018       if (indphi.gt.0) then
1019         nvar=nphi
1020       else if (indback.gt.0) then
1021         nvar=nphi+ntheta
1022       else
1023         nvar=nvar+2*nside
1024       endif
1025 cd    write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1026       return
1027       end
1028 c----------------------------------------------------------------------------
1029       subroutine gen_dist_constr
1030 C Generate CA distance constraints.
1031       implicit real*8 (a-h,o-z)
1032       include 'DIMENSIONS'
1033       include 'COMMON.IOUNITS'
1034       include 'COMMON.GEO'
1035       include 'COMMON.VAR'
1036       include 'COMMON.INTERACT'
1037       include 'COMMON.LOCAL'
1038       include 'COMMON.NAMES'
1039       include 'COMMON.CHAIN'
1040       include 'COMMON.FFIELD'
1041       include 'COMMON.SBRIDGE'
1042       include 'COMMON.HEADER'
1043       include 'COMMON.CONTROL'
1044 c      include 'COMMON.DBASE'
1045 c      include 'COMMON.THREAD'
1046       include 'COMMON.TIME1'
1047       dimension itype_pdb(maxres)
1048       common /pizda/ itype_pdb
1049       character*2 iden
1050 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1051 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1052 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1053 cd     & ' nsup',nsup
1054       do i=nstart_sup,nstart_sup+nsup-1
1055 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1056 cd     &    ' seq_pdb', restyp(itype_pdb(i))
1057         do j=i+2,nstart_sup+nsup-1
1058           nhpb=nhpb+1
1059           ihpb(nhpb)=i+nstart_seq-nstart_sup
1060           jhpb(nhpb)=j+nstart_seq-nstart_sup
1061           forcon(nhpb)=weidis
1062           dhpb(nhpb)=dist(i,j)
1063         enddo
1064       enddo 
1065 cd      write (iout,'(a)') 'Distance constraints:' 
1066 cd      do i=nss+1,nhpb
1067 cd        ii=ihpb(i)
1068 cd        jj=jhpb(i)
1069 cd        iden='CA'
1070 cd        if (ii.gt.nres) then
1071 cd          iden='SC'
1072 cd          ii=ii-nres
1073 cd          jj=jj-nres
1074 cd        endif
1075 cd        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)') 
1076 cd     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1077 cd     &  dhpb(i),forcon(i)
1078 cd      enddo
1079       return
1080       end
1081 c----------------------------------------------------------------------------
1082       subroutine csaread
1083       implicit real*8 (a-h,o-z)
1084       include 'DIMENSIONS'
1085       include 'COMMON.IOUNITS'
1086       include 'COMMON.GEO'
1087       include 'COMMON.CSA'
1088       include 'COMMON.BANK'
1089       include 'COMMON.CONTROL'
1090       character*80 ucase
1091       character*620 mcmcard
1092       call card_concat(mcmcard)
1093
1094       call readi(mcmcard,'NCONF',nconf,50)
1095       call readi(mcmcard,'NADD',nadd,0)
1096       call readi(mcmcard,'JSTART',jstart,1)
1097       call readi(mcmcard,'JEND',jend,1)
1098       call readi(mcmcard,'NSTMAX',nstmax,500000)
1099       call readi(mcmcard,'N0',n0,1)
1100       call readi(mcmcard,'N1',n1,6)
1101       call readi(mcmcard,'N2',n2,4)
1102       call readi(mcmcard,'N3',n3,0)
1103       call readi(mcmcard,'N4',n4,0)
1104       call readi(mcmcard,'N5',n5,0)
1105       call readi(mcmcard,'N6',n6,10)
1106       call readi(mcmcard,'N7',n7,0)
1107       call readi(mcmcard,'N8',n8,0)
1108       call readi(mcmcard,'N9',n9,0)
1109       call readi(mcmcard,'N14',n14,0)
1110       call readi(mcmcard,'N15',n15,0)
1111       call readi(mcmcard,'N16',n16,0)
1112       call readi(mcmcard,'N17',n17,0)
1113       call readi(mcmcard,'N18',n18,0)
1114
1115       vdisulf=(index(mcmcard,'DYNSS').gt.0)
1116
1117       call readi(mcmcard,'NDIFF',ndiff,2)
1118       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1119       call readi(mcmcard,'IS1',is1,1)
1120       call readi(mcmcard,'IS2',is2,8)
1121       call readi(mcmcard,'NRAN0',nran0,4)
1122       call readi(mcmcard,'NRAN1',nran1,2)
1123       call readi(mcmcard,'IRR',irr,1)
1124       call readi(mcmcard,'NSEED',nseed,20)
1125       call readi(mcmcard,'NTOTAL',ntotal,10000)
1126       call reada(mcmcard,'CUT1',cut1,2.0d0)
1127       call reada(mcmcard,'CUT2',cut2,5.0d0)
1128       call reada(mcmcard,'ESTOP',estop,-300000.0d0)
1129       call readi(mcmcard,'ICMAX',icmax,1)
1130       call readi(mcmcard,'NBANKM',nbankm,400)
1131       call readi(mcmcard,'IUCUT',iucut,2)
1132       call readi(mcmcard,'IRESTART',irestart,0)
1133 c!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
1134       ntbankm=0
1135 c!bankt
1136       call reada(mcmcard,'DELE',dele,20.0d0)
1137       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1138       call readi(mcmcard,'IREF',iref,0)
1139       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1140       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1141       call readi(mcmcard,'NCONF_IN',nconf_in,0)
1142       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1143       write (iout,*) "NCONF_IN",nconf_in
1144       tm_score=(index(mcmcard,'TMSCORE').gt.0)
1145       if (tm_score) write (iout,*) "Using TM_Score instead of DIFF",
1146      &      " for torsional angles" 
1147       return
1148       end
1149
1150       subroutine read_minim
1151       implicit real*8 (a-h,o-z)
1152       include 'DIMENSIONS'
1153       include 'COMMON.MINIM'
1154       include 'COMMON.IOUNITS'
1155       character*80 ucase
1156       character*320 minimcard
1157       call card_concat(minimcard)
1158       call readi(minimcard,'MAXMIN',maxmin,2000)
1159       call readi(minimcard,'MAXFUN',maxfun,5000)
1160       call readi(minimcard,'MINMIN',minmin,maxmin)
1161       call readi(minimcard,'MINFUN',minfun,maxmin)
1162       call reada(minimcard,'TOLF',tolf,1.0D-2)
1163       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1164       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1165       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1166       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1167       write (iout,'(/80(1h*)/20x,a/80(1h*))') 
1168      &         'Options in energy minimization:'
1169       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1170      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1171      & 'MinMin:',MinMin,' MinFun:',MinFun,
1172      & ' TolF:',TolF,' RTolF:',RTolF
1173       return
1174       end
1175 c----------------------------------------------------------------------------
1176       subroutine read_angles(kanal,*)
1177       implicit real*8 (a-h,o-z)
1178       include 'DIMENSIONS'
1179       include 'COMMON.GEO'
1180       include 'COMMON.VAR'
1181       include 'COMMON.CHAIN'
1182       include 'COMMON.IOUNITS'
1183       include 'COMMON.CONTROL'
1184 c Read angles from input 
1185 c
1186        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1187        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1188        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1189        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1190
1191        do i=1,nres
1192 c 9/7/01 avoid 180 deg valence angle
1193         if (theta(i).gt.179.99d0) theta(i)=179.99d0
1194 c
1195         theta(i)=deg2rad*theta(i)
1196         phi(i)=deg2rad*phi(i)
1197         alph(i)=deg2rad*alph(i)
1198         omeg(i)=deg2rad*omeg(i)
1199        enddo
1200       return
1201    10 return1
1202       end
1203 c----------------------------------------------------------------------------
1204       subroutine reada(rekord,lancuch,wartosc,default)
1205       implicit none
1206       character*(*) rekord,lancuch
1207       double precision wartosc,default
1208       integer ilen,iread
1209       external ilen
1210       iread=index(rekord,lancuch)
1211       if (iread.eq.0) then
1212         wartosc=default 
1213         return
1214       endif   
1215       iread=iread+ilen(lancuch)+1
1216       read (rekord(iread:),*,err=10,end=10) wartosc
1217       return
1218   10  wartosc=default
1219       return
1220       end
1221 c----------------------------------------------------------------------------
1222       subroutine readi(rekord,lancuch,wartosc,default)
1223       implicit none
1224       character*(*) rekord,lancuch
1225       integer wartosc,default
1226       integer ilen,iread
1227       external ilen
1228       iread=index(rekord,lancuch)
1229       if (iread.eq.0) then
1230         wartosc=default 
1231         return
1232       endif   
1233       iread=iread+ilen(lancuch)+1
1234       read (rekord(iread:),*,err=10,end=10) wartosc
1235       return
1236   10  wartosc=default
1237       return
1238       end
1239 c----------------------------------------------------------------------------
1240       subroutine multreadi(rekord,lancuch,tablica,dim,default)
1241       implicit none
1242       integer dim,i
1243       integer tablica(dim),default
1244       character*(*) rekord,lancuch
1245       character*80 aux
1246       integer ilen,iread
1247       external ilen
1248       do i=1,dim
1249         tablica(i)=default
1250       enddo
1251       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1252       if (iread.eq.0) return
1253       iread=iread+ilen(lancuch)+1
1254       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1255    10 return
1256       end
1257 c----------------------------------------------------------------------------
1258       subroutine multreada(rekord,lancuch,tablica,dim,default)
1259       implicit none
1260       integer dim,i
1261       double precision tablica(dim),default
1262       character*(*) rekord,lancuch
1263       character*80 aux
1264       integer ilen,iread
1265       external ilen
1266       do i=1,dim
1267         tablica(i)=default
1268       enddo
1269       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1270       if (iread.eq.0) return
1271       iread=iread+ilen(lancuch)+1
1272       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1273    10 return
1274       end
1275 c----------------------------------------------------------------------------
1276       subroutine openunits
1277       implicit real*8 (a-h,o-z)
1278       include 'DIMENSIONS'    
1279 #ifdef MPI
1280       include 'mpif.h'
1281       character*16 form,nodename
1282       integer nodelen
1283 #endif
1284       include 'COMMON.SETUP'
1285       include 'COMMON.IOUNITS'
1286 c      include 'COMMON.MD'
1287       include 'COMMON.CONTROL'
1288       integer lenpre,lenpot,ilen,lentmp
1289       external ilen
1290       character*3 out1file_text,ucase
1291       character*3 ll
1292       external ucase
1293 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1294       call getenv_loc("PREFIX",prefix)
1295       pref_orig = prefix
1296       call getenv_loc("POT",pot)
1297       call getenv_loc("DIRTMP",tmpdir)
1298       call getenv_loc("CURDIR",curdir)
1299       call getenv_loc("OUT1FILE",out1file_text)
1300 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1301       out1file_text=ucase(out1file_text)
1302       if (out1file_text(1:1).eq."Y") then
1303         out1file=.true.
1304       else 
1305         out1file=fg_rank.gt.0
1306       endif
1307       lenpre=ilen(prefix)
1308       lenpot=ilen(pot)
1309       lentmp=ilen(tmpdir)
1310       if (lentmp.gt.0) then
1311           write (*,'(80(1h!))')
1312           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
1313           write (*,'(80(1h!))')
1314           write (*,*)"All output files will be on node /tmp directory." 
1315 #ifdef MPI
1316         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1317         if (me.eq.king) then
1318           write (*,*) "The master node is ",nodename
1319         else if (fg_rank.eq.0) then
1320           write (*,*) "I am the CG slave node ",nodename
1321         else 
1322           write (*,*) "I am the FG slave node ",nodename
1323         endif
1324 #endif
1325         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1326         lenpre = lentmp+lenpre+1
1327       endif
1328       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1329 C Get the names and open the input files
1330 #if defined(WINIFL) || defined(WINPGI)
1331       open(1,file=pref_orig(:ilen(pref_orig))//
1332      &  '.inp',status='old',readonly,shared)
1333        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1334 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1335 C Get parameter filenames and open the parameter files.
1336       call getenv_loc('BONDPAR',bondname)
1337       open (ibond,file=bondname,status='old',readonly,shared)
1338       call getenv_loc('THETPAR',thetname)
1339       open (ithep,file=thetname,status='old',readonly,shared)
1340 #ifndef CRYST_THETA
1341       call getenv_loc('THETPARPDB',thetname_pdb)
1342       open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
1343 #endif
1344       call getenv_loc('ROTPAR',rotname)
1345       open (irotam,file=rotname,status='old',readonly,shared)
1346 #ifndef CRYST_SC
1347       call getenv_loc('ROTPARPDB',rotname_pdb)
1348       open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
1349 #endif
1350       call getenv_loc('TORPAR',torname)
1351       open (itorp,file=torname,status='old',readonly,shared)
1352       call getenv_loc('TORDPAR',tordname)
1353       open (itordp,file=tordname,status='old',readonly,shared)
1354       call getenv_loc('FOURIER',fouriername)
1355       open (ifourier,file=fouriername,status='old',readonly,shared)
1356       call getenv_loc('ELEPAR',elename)
1357       open (ielep,file=elename,status='old',readonly,shared)
1358       call getenv_loc('SIDEPAR',sidename)
1359       open (isidep,file=sidename,status='old',readonly,shared)
1360 #elif (defined CRAY) || (defined AIX)
1361       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1362      &  action='read')
1363 c      print *,"Processor",myrank," opened file 1" 
1364       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1365 c      print *,"Processor",myrank," opened file 9" 
1366 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1367 C Get parameter filenames and open the parameter files.
1368       call getenv_loc('BONDPAR',bondname)
1369       open (ibond,file=bondname,status='old',action='read')
1370 c      print *,"Processor",myrank," opened file IBOND" 
1371       call getenv_loc('THETPAR',thetname)
1372       open (ithep,file=thetname,status='old',action='read')
1373 c      print *,"Processor",myrank," opened file ITHEP" 
1374 #ifndef CRYST_THETA
1375       call getenv_loc('THETPARPDB',thetname_pdb)
1376       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1377 #endif
1378       call getenv_loc('ROTPAR',rotname)
1379       open (irotam,file=rotname,status='old',action='read')
1380 c      print *,"Processor",myrank," opened file IROTAM" 
1381 #ifndef CRYST_SC
1382       call getenv_loc('ROTPARPDB',rotname_pdb)
1383       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
1384 #endif
1385       call getenv_loc('TORPAR',torname)
1386       open (itorp,file=torname,status='old',action='read')
1387 c      print *,"Processor",myrank," opened file ITORP" 
1388       call getenv_loc('TORDPAR',tordname)
1389       open (itordp,file=tordname,status='old',action='read')
1390 c      print *,"Processor",myrank," opened file ITORDP" 
1391       call getenv_loc('SCCORPAR',sccorname)
1392       open (isccor,file=sccorname,status='old',action='read')
1393 c      print *,"Processor",myrank," opened file ISCCOR" 
1394       call getenv_loc('FOURIER',fouriername)
1395       open (ifourier,file=fouriername,status='old',action='read')
1396 c      print *,"Processor",myrank," opened file IFOURIER" 
1397       call getenv_loc('ELEPAR',elename)
1398       open (ielep,file=elename,status='old',action='read')
1399 c      print *,"Processor",myrank," opened file IELEP" 
1400       call getenv_loc('SIDEPAR',sidename)
1401       open (isidep,file=sidename,status='old',action='read')
1402 c      print *,"Processor",myrank," opened file ISIDEP" 
1403 c      print *,"Processor",myrank," opened parameter files" 
1404 #elif (defined G77)
1405       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1406       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1407 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1408 C Get parameter filenames and open the parameter files.
1409       call getenv_loc('BONDPAR',bondname)
1410       open (ibond,file=bondname,status='old')
1411       call getenv_loc('THETPAR',thetname)
1412       open (ithep,file=thetname,status='old')
1413 #ifndef CRYST_THETA
1414       call getenv_loc('THETPARPDB',thetname_pdb)
1415       open (ithep_pdb,file=thetname_pdb,status='old')
1416 #endif
1417       call getenv_loc('ROTPAR',rotname)
1418       open (irotam,file=rotname,status='old')
1419 #ifndef CRYST_SC
1420       call getenv_loc('ROTPARPDB',rotname_pdb)
1421       open (irotam_pdb,file=rotname_pdb,status='old')
1422 #endif
1423       call getenv_loc('TORPAR',torname)
1424       open (itorp,file=torname,status='old')
1425       call getenv_loc('TORDPAR',tordname)
1426       open (itordp,file=tordname,status='old')
1427       call getenv_loc('SCCORPAR',sccorname)
1428       open (isccor,file=sccorname,status='old')
1429       call getenv_loc('FOURIER',fouriername)
1430       open (ifourier,file=fouriername,status='old')
1431       call getenv_loc('ELEPAR',elename)
1432       open (ielep,file=elename,status='old')
1433       call getenv_loc('SIDEPAR',sidename)
1434       open (isidep,file=sidename,status='old')
1435 #else
1436       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1437      &  readonly)
1438        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1439 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1440 C Get parameter filenames and open the parameter files.
1441       call getenv_loc('BONDPAR',bondname)
1442       open (ibond,file=bondname,status='old',readonly)
1443       call getenv_loc('THETPAR',thetname)
1444       open (ithep,file=thetname,status='old',readonly)
1445 #ifndef CRYST_THETA
1446       call getenv_loc('THETPARPDB',thetname_pdb)
1447       print *,"thetname_pdb ",thetname_pdb
1448       open (ithep_pdb,file=thetname_pdb,status='old',readonly)
1449       print *,ithep_pdb," opened"
1450 #endif
1451       call getenv_loc('ROTPAR',rotname)
1452       open (irotam,file=rotname,status='old',readonly)
1453 #ifndef CRYST_SC
1454       call getenv_loc('ROTPARPDB',rotname_pdb)
1455       open (irotam_pdb,file=rotname_pdb,status='old',readonly)
1456 #endif
1457       call getenv_loc('TORPAR',torname)
1458       open (itorp,file=torname,status='old',readonly)
1459       call getenv_loc('TORDPAR',tordname)
1460       open (itordp,file=tordname,status='old',readonly)
1461       call getenv_loc('SCCORPAR',sccorname)
1462       open (isccor,file=sccorname,status='old',readonly)
1463       call getenv_loc('FOURIER',fouriername)
1464       open (ifourier,file=fouriername,status='old',readonly)
1465       call getenv_loc('ELEPAR',elename)
1466       open (ielep,file=elename,status='old',readonly)
1467       call getenv_loc('SIDEPAR',sidename)
1468       open (isidep,file=sidename,status='old',readonly)
1469 #endif
1470 #ifndef OLDSCP
1471 C
1472 C 8/9/01 In the newest version SCp interaction constants are read from a file
1473 C Use -DOLDSCP to use hard-coded constants instead.
1474 C
1475       call getenv_loc('SCPPAR',scpname)
1476 #if defined(WINIFL) || defined(WINPGI)
1477       open (iscpp,file=scpname,status='old',readonly,shared)
1478 #elif (defined CRAY)  || (defined AIX)
1479       open (iscpp,file=scpname,status='old',action='read')
1480 #elif (defined G77)
1481       open (iscpp,file=scpname,status='old')
1482 #else
1483       open (iscpp,file=scpname,status='old',readonly)
1484 #endif
1485 #endif
1486       call getenv_loc('PATTERN',patname)
1487 #if defined(WINIFL) || defined(WINPGI)
1488       open (icbase,file=patname,status='old',readonly,shared)
1489 #elif (defined CRAY)  || (defined AIX)
1490       open (icbase,file=patname,status='old',action='read')
1491 #elif (defined G77)
1492       open (icbase,file=patname,status='old')
1493 #else
1494       open (icbase,file=patname,status='old',readonly)
1495 #endif
1496 #ifdef MPI
1497 C Open output file only for CG processes
1498 c      print *,"Processor",myrank," fg_rank",fg_rank
1499       if (fg_rank.eq.0) then
1500
1501       if (nodes.eq.1) then
1502         npos=3
1503       else
1504         npos = dlog10(dfloat(nodes-1))+1
1505       endif
1506       if (npos.lt.3) npos=3
1507       write (liczba,'(i1)') npos
1508       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
1509      &  //')'
1510       write (liczba,form) me
1511       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
1512      &  liczba(:ilen(liczba))
1513       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1514      &  //'.int'
1515       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1516      &  //'.pdb'
1517       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
1518      &  liczba(:ilen(liczba))//'.mol2'
1519       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
1520      &  liczba(:ilen(liczba))//'.stat'
1521       if (lentmp.gt.0)
1522      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
1523      &      //liczba(:ilen(liczba))//'.stat')
1524       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
1525      &  //'.rst'
1526 c      if(usampl) then
1527 c          qname=prefix(:lenpre)//'_'//pot(:lenpot)//
1528 c     & liczba(:ilen(liczba))//'.const'
1529 c      endif 
1530
1531       endif
1532 #else
1533       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
1534       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
1535       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
1536       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
1537       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
1538       if (lentmp.gt.0)
1539      & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)//
1540      &    '.stat')
1541       rest2name=prefix(:ilen(prefix))//'.rst'
1542 c      if(usampl) then 
1543 c         qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
1544 c      endif 
1545 #endif
1546 #if defined(AIX) || defined(PGI)
1547       if (me.eq.king .or. .not. out1file) 
1548      &   open(iout,file=outname,status='unknown')
1549 #ifdef DEBUG
1550       if (fg_rank.gt.0) then
1551         write (liczba,'(i3.3)') myrank/nfgtasks
1552         write (ll,'(bz,i3.3)') fg_rank
1553         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1554      &   status='unknown')
1555       endif
1556 #endif
1557       if(me.eq.king) then
1558        open(igeom,file=intname,status='unknown',position='append')
1559        open(ipdb,file=pdbname,status='unknown')
1560        open(imol2,file=mol2name,status='unknown')
1561        open(istat,file=statname,status='unknown',position='append')
1562       else
1563 c1out       open(iout,file=outname,status='unknown')
1564       endif
1565 #else
1566       if (me.eq.king .or. .not.out1file)
1567      &    open(iout,file=outname,status='unknown')
1568 #ifdef DEBUG
1569       if (fg_rank.gt.0) then
1570         write (liczba,'(i3.3)') myrank/nfgtasks
1571         write (ll,'(bz,i3.3)') fg_rank
1572         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1573      &   status='unknown')
1574       endif
1575 #endif
1576       if(me.eq.king) then
1577        open(igeom,file=intname,status='unknown',access='append')
1578        open(ipdb,file=pdbname,status='unknown')
1579        open(imol2,file=mol2name,status='unknown')
1580        open(istat,file=statname,status='unknown',access='append')
1581       else
1582 c1out       open(iout,file=outname,status='unknown')
1583       endif
1584 #endif
1585       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
1586       csa_seed=prefix(:lenpre)//'.CSA.seed'
1587       csa_history=prefix(:lenpre)//'.CSA.history'
1588       csa_bank=prefix(:lenpre)//'.CSA.bank'
1589       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
1590       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
1591       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
1592 c!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
1593       csa_int=prefix(:lenpre)//'.int'
1594       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
1595       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
1596       csa_in=prefix(:lenpre)//'.CSA.in'
1597 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
1598 C Write file names
1599       if (me.eq.king)then
1600       write (iout,'(80(1h-))')
1601       write (iout,'(30x,a)') "FILE ASSIGNMENT"
1602       write (iout,'(80(1h-))')
1603       write (iout,*) "Input file                      : ",
1604      &  pref_orig(:ilen(pref_orig))//'.inp'
1605       write (iout,*) "Output file                     : ",
1606      &  outname(:ilen(outname))
1607       write (iout,*)
1608       write (iout,*) "Sidechain potential file        : ",
1609      &  sidename(:ilen(sidename))
1610 #ifndef OLDSCP
1611       write (iout,*) "SCp potential file              : ",
1612      &  scpname(:ilen(scpname))
1613 #endif
1614       write (iout,*) "Electrostatic potential file    : ",
1615      &  elename(:ilen(elename))
1616       write (iout,*) "Cumulant coefficient file       : ",
1617      &  fouriername(:ilen(fouriername))
1618       write (iout,*) "Torsional parameter file        : ",
1619      &  torname(:ilen(torname))
1620       write (iout,*) "Double torsional parameter file : ",
1621      &  tordname(:ilen(tordname))
1622       write (iout,*) "SCCOR parameter file : ",
1623      &  sccorname(:ilen(sccorname))
1624       write (iout,*) "Bond & inertia constant file    : ",
1625      &  bondname(:ilen(bondname))
1626       write (iout,*) "Bending parameter file          : ",
1627      &  thetname(:ilen(thetname))
1628       write (iout,*) "Rotamer parameter file          : ",
1629      &  rotname(:ilen(rotname))
1630       write (iout,*) "Threading database              : ",
1631      &  patname(:ilen(patname))
1632       if (lentmp.ne.0) 
1633      &write (iout,*)" DIRTMP                          : ",
1634      &  tmpdir(:lentmp)
1635       write (iout,'(80(1h-))')
1636       endif
1637       return
1638       end
1639 c----------------------------------------------------------------------------
1640       subroutine card_concat(card)
1641       implicit real*8 (a-h,o-z)
1642       include 'DIMENSIONS'
1643       include 'COMMON.IOUNITS'
1644       character*(*) card
1645       character*80 karta,ucase
1646       external ilen
1647       read (inp,'(a)') karta
1648       karta=ucase(karta)
1649       card=' '
1650       do while (karta(80:80).eq.'&')
1651         card=card(:ilen(card)+1)//karta(:79)
1652         read (inp,'(a)') karta
1653         karta=ucase(karta)
1654       enddo
1655       card=card(:ilen(card)+1)//karta
1656       return
1657       end
1658
1659       subroutine read_dist_constr
1660       implicit real*8 (a-h,o-z)
1661       include 'DIMENSIONS'
1662 #ifdef MPI
1663       include 'mpif.h'
1664 #endif
1665       include 'COMMON.SETUP'
1666       include 'COMMON.CONTROL'
1667       include 'COMMON.CHAIN'
1668       include 'COMMON.IOUNITS'
1669       include 'COMMON.SBRIDGE'
1670       integer ifrag_(2,100),ipair_(2,100)
1671       double precision wfrag_(100),wpair_(100)
1672       character*500 controlcard
1673 c      write (iout,*) "Calling read_dist_constr"
1674 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
1675 c      call flush(iout)
1676       call card_concat(controlcard)
1677       call readi(controlcard,"NFRAG",nfrag_,0)
1678       call readi(controlcard,"NPAIR",npair_,0)
1679       call readi(controlcard,"NDIST",ndist_,0)
1680       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
1681       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
1682       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
1683       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
1684       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
1685 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
1686 c      write (iout,*) "IFRAG"
1687 c      do i=1,nfrag_
1688 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1689 c      enddo
1690 c      write (iout,*) "IPAIR"
1691 c      do i=1,npair_
1692 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
1693 c      enddo
1694       if (.not.refstr .and. nfrag.gt.0) then
1695         write (iout,*) 
1696      &  "ERROR: no reference structure to compute distance restraints"
1697         write (iout,*)
1698      &  "Restraints must be specified explicitly (NDIST=number)"
1699         stop 
1700       endif
1701       if (nfrag.lt.2 .and. npair.gt.0) then 
1702         write (iout,*) "ERROR: Less than 2 fragments specified",
1703      &   " but distance restraints between pairs requested"
1704         stop 
1705       endif 
1706       call flush(iout)
1707       do i=1,nfrag_
1708         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
1709         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
1710      &    ifrag_(2,i)=nstart_sup+nsup-1
1711 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1712         call flush(iout)
1713         if (wfrag_(i).gt.0.0d0) then
1714         do j=ifrag_(1,i),ifrag_(2,i)-1
1715           do k=j+1,ifrag_(2,i)
1716 c            write (iout,*) "j",j," k",k
1717             ddjk=dist(j,k)
1718             if (constr_dist.eq.1) then
1719             nhpb=nhpb+1
1720             ihpb(nhpb)=j
1721             jhpb(nhpb)=k
1722               dhpb(nhpb)=ddjk
1723             forcon(nhpb)=wfrag_(i) 
1724             else if (constr_dist.eq.2) then
1725               if (ddjk.le.dist_cut) then
1726                 nhpb=nhpb+1
1727                 ihpb(nhpb)=j
1728                 jhpb(nhpb)=k
1729                 dhpb(nhpb)=ddjk
1730                 forcon(nhpb)=wfrag_(i) 
1731               endif
1732             else
1733               nhpb=nhpb+1
1734               ihpb(nhpb)=j
1735               jhpb(nhpb)=k
1736               dhpb(nhpb)=ddjk
1737               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1738             endif
1739 #ifdef MPI
1740             if (.not.out1file .or. me.eq.king) 
1741      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1742      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1743 #else
1744             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1745      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1746 #endif
1747           enddo
1748         enddo
1749         endif
1750       enddo
1751       do i=1,npair_
1752         if (wpair_(i).gt.0.0d0) then
1753         ii = ipair_(1,i)
1754         jj = ipair_(2,i)
1755         if (ii.gt.jj) then
1756           itemp=ii
1757           ii=jj
1758           jj=itemp
1759         endif
1760         do j=ifrag_(1,ii),ifrag_(2,ii)
1761           do k=ifrag_(1,jj),ifrag_(2,jj)
1762             nhpb=nhpb+1
1763             ihpb(nhpb)=j
1764             jhpb(nhpb)=k
1765             forcon(nhpb)=wpair_(i)
1766             dhpb(nhpb)=dist(j,k)
1767 #ifdef MPI
1768             if (.not.out1file .or. me.eq.king)
1769      &      write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1770      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1771 #else
1772             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1773      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1774 #endif
1775           enddo
1776         enddo
1777         endif
1778       enddo 
1779       do i=1,ndist_
1780         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
1781      &     ibecarb(i),forcon(nhpb+1)
1782         if (forcon(nhpb+1).gt.0.0d0) then
1783           nhpb=nhpb+1
1784           if (ibecarb(i).gt.0) then
1785             ihpb(i)=ihpb(i)+nres
1786             jhpb(i)=jhpb(i)+nres
1787           endif
1788           if (dhpb(nhpb).eq.0.0d0) 
1789      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1790         endif
1791       enddo
1792 #ifdef MPI
1793       if (.not.out1file .or. me.eq.king) then
1794 #endif
1795       do i=1,nhpb
1796           write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
1797      &     i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
1798       enddo
1799       call flush(iout)
1800 #ifdef MPI
1801       endif
1802 #endif
1803       return
1804       end
1805 c-------------------------------------------------------------------------------
1806 #ifdef WINIFL
1807       subroutine flush(iu)
1808       return
1809       end
1810 #endif
1811 #ifdef AIX
1812       subroutine flush(iu)
1813       call flush_(iu)
1814       return
1815       end
1816 #endif
1817 c------------------------------------------------------------------------------
1818       subroutine copy_to_tmp(source)
1819       include "DIMENSIONS"
1820       include "COMMON.IOUNITS"
1821       character*(*) source
1822       character* 256 tmpfile
1823       integer ilen
1824       external ilen
1825       logical ex
1826       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
1827       inquire(file=tmpfile,exist=ex)
1828       if (ex) then
1829         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
1830      &   " to temporary directory..."
1831         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
1832         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
1833       endif
1834       return
1835       end
1836 c------------------------------------------------------------------------------
1837       subroutine move_from_tmp(source)
1838       include "DIMENSIONS"
1839       include "COMMON.IOUNITS"
1840       character*(*) source
1841       integer ilen
1842       external ilen
1843       write (*,*) "Moving ",source(:ilen(source)),
1844      & " from temporary directory to working directory"
1845       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
1846       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
1847       return
1848       end
1849 c------------------------------------------------------------------------------
1850       subroutine random_init(seed)
1851 C
1852 C Initialize random number generator
1853 C
1854       implicit real*8 (a-h,o-z)
1855       include 'DIMENSIONS'
1856 #ifdef AMD64
1857       integer*8 iseedi8
1858 #endif
1859 #ifdef MPI
1860       include 'mpif.h'
1861       logical OKRandom, prng_restart
1862       real*8  r1
1863       integer iseed_array(4)
1864 #endif
1865       include 'COMMON.IOUNITS'
1866       include 'COMMON.TIME1'
1867 c      include 'COMMON.THREAD'
1868       include 'COMMON.SBRIDGE'
1869       include 'COMMON.CONTROL'
1870       include 'COMMON.MCM'
1871 c      include 'COMMON.MAP'
1872       include 'COMMON.HEADER'
1873 c      include 'COMMON.CSA'
1874       include 'COMMON.CHAIN'
1875 c      include 'COMMON.MUCA'
1876 c      include 'COMMON.MD'
1877       include 'COMMON.FFIELD'
1878       include 'COMMON.SETUP'
1879       iseed=-dint(dabs(seed))
1880       if (iseed.eq.0) then
1881         write (iout,'(/80(1h*)/20x,a/80(1h*))') 
1882      &    'Random seed undefined. The program will stop.'
1883         write (*,'(/80(1h*)/20x,a/80(1h*))') 
1884      &    'Random seed undefined. The program will stop.'
1885 #ifdef MPI
1886         call mpi_finalize(mpi_comm_world,ierr)
1887 #endif
1888         stop 'Bad random seed.'
1889       endif
1890 #ifdef MPI
1891       if (fg_rank.eq.0) then
1892       seed=seed*(me+1)+1
1893 #ifdef AMD64
1894       iseedi8=dint(seed)
1895       if(me.eq.king .or. .not. out1file)
1896      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1897       write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1898       OKRandom = prng_restart(me,iseedi8)
1899 #else
1900       do i=1,4
1901        tmp=65536.0d0**(4-i)
1902        iseed_array(i) = dint(seed/tmp)
1903        seed=seed-iseed_array(i)*tmp
1904       enddo
1905       if(me.eq.king .or. .not. out1file)
1906      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
1907      &                 (iseed_array(i),i=1,4)
1908       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
1909      &                 (iseed_array(i),i=1,4)
1910       OKRandom = prng_restart(me,iseed_array)
1911 #endif
1912       if (OKRandom) then
1913         r1=ran_number(0.0D0,1.0D0)
1914         if(me.eq.king .or. .not. out1file)
1915      &   write (iout,*) 'ran_num',r1
1916         if (r1.lt.0.0d0) OKRandom=.false.
1917       endif
1918       if (.not.OKRandom) then
1919         write (iout,*) 'PRNG IS NOT WORKING!!!'
1920         print *,'PRNG IS NOT WORKING!!!'
1921         if (me.eq.0) then 
1922          call flush(iout)
1923          call mpi_abort(mpi_comm_world,error_msg,ierr)
1924          stop
1925         else
1926          write (iout,*) 'too many processors for parallel prng'
1927          write (*,*) 'too many processors for parallel prng'
1928          call flush(iout)
1929          stop
1930         endif
1931       endif
1932       endif
1933 #else
1934       call vrndst(iseed)
1935       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
1936 #endif
1937       return
1938       end