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