2 implicit real*8 (a-h,o-z)
8 include 'COMMON.CONTROL'
9 include 'COMMON.SBRIDGE'
10 include 'COMMON.IOUNITS'
11 include 'COMMON.SPLITELE'
13 C Read force-field parameters except weights
15 C Read job setup parameters
17 C Read control parameters for energy minimzation if required
18 if (minim) call read_minim
19 C Read MCM control parameters if required
20 if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
21 C Read MD control parameters if reqjuired
22 if (modecalc.eq.12) call read_MDpar
23 C Read MREMD control parameters if required
24 if (modecalc.eq.14) then
28 C Read MUCA control parameters if required
29 if (lmuca) call read_muca
30 C Read CSA control parameters if required (from fort.40 if exists
31 C otherwise from general input file)
32 if (modecalc.eq.8) then
33 inquire (file="fort.40",exist=file_exist)
34 if (.not.file_exist) call csaread
36 cfmc if (modecalc.eq.10) call mcmfread
37 C Read molecule information, molecule geometry, energy-term weights, and
38 C restraints if requested
40 C Print restraint information
42 if (.not. out1file .or. me.eq.king) then
45 &write (iout,'(a,i5,a)') "The following",nhpb-nss,
46 & " distance constraints have been imposed"
48 write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
53 c print *,"Processor",myrank," leaves READRTNS"
56 C-------------------------------------------------------------------------------
57 subroutine read_control
61 implicit real*8 (a-h,o-z)
65 logical OKRandom, prng_restart
68 include 'COMMON.IOUNITS'
69 include 'COMMON.TIME1'
70 include 'COMMON.THREAD'
71 include 'COMMON.SBRIDGE'
72 include 'COMMON.CONTROL'
75 include 'COMMON.HEADER'
77 include 'COMMON.CHAIN'
80 include 'COMMON.FFIELD'
81 include 'COMMON.INTERACT'
82 include 'COMMON.SETUP'
83 include 'COMMON.SPLITELE'
84 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
85 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
87 character*320 controlcard
92 read (INP,'(a)') titel
93 call card_concat(controlcard)
94 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
95 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
96 call reada(controlcard,'SEED',seed,0.0D0)
97 call random_init(seed)
98 C Set up the time limit (caution! The time must be input in minutes!)
99 read_cart=index(controlcard,'READ_CART').gt.0
100 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
101 C this variable with_theta_constr is the variable which allow to read and execute the
102 C constrains on theta angles WITH_THETA_CONSTR is the keyword
103 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
104 write (iout,*) "with_theta_constr ",with_theta_constr
105 call readi(controlcard,'SYM',symetr,1)
106 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
107 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
108 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
109 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
110 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
111 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
112 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
113 call reada(controlcard,'DRMS',drms,0.1D0)
114 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
115 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
116 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
117 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
118 write (iout,'(a,f10.1)')'DRMS = ',drms
119 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
120 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
122 call readi(controlcard,'NZ_START',nz_start,0)
123 call readi(controlcard,'NZ_END',nz_end,0)
124 call readi(controlcard,'IZ_SC',iz_sc,0)
126 safety = 60.0d0*safety
129 call reada(controlcard,"T_BATH",t_bath,300.0d0)
130 minim=(index(controlcard,'MINIMIZE').gt.0)
131 dccart=(index(controlcard,'CART').gt.0)
132 overlapsc=(index(controlcard,'OVERLAP').gt.0)
133 overlapsc=.not.overlapsc
134 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
135 searchsc=.not.searchsc
136 sideadd=(index(controlcard,'SIDEADD').gt.0)
137 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
138 outpdb=(index(controlcard,'PDBOUT').gt.0)
139 outmol2=(index(controlcard,'MOL2OUT').gt.0)
140 pdbref=(index(controlcard,'PDBREF').gt.0)
141 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
142 indpdb=index(controlcard,'PDBSTART')
143 extconf=(index(controlcard,'EXTCONF').gt.0)
144 AFMlog=(index(controlcard,'AFM'))
145 selfguide=(index(controlcard,'SELFGUIDE'))
146 print *,'AFMlog',AFMlog,selfguide,"KUPA"
147 call readi(controlcard,'IPRINT',iprint,0)
148 call readi(controlcard,'MAXGEN',maxgen,10000)
149 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
150 call readi(controlcard,"KDIAG",kdiag,0)
151 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
152 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
153 & write (iout,*) "RESCALE_MODE",rescale_mode
154 split_ene=index(controlcard,'SPLIT_ENE').gt.0
155 if (index(controlcard,'REGULAR').gt.0.0D0) then
156 call reada(controlcard,'WEIDIS',weidis,0.1D0)
160 if (index(controlcard,'CHECKGRAD').gt.0) then
162 if (index(controlcard,'CART').gt.0) then
164 elseif (index(controlcard,'CARINT').gt.0) then
169 elseif (index(controlcard,'THREAD').gt.0) then
171 call readi(controlcard,'THREAD',nthread,0)
172 if (nthread.gt.0) then
173 call reada(controlcard,'WEIDIS',weidis,0.1D0)
176 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
177 stop 'Error termination in Read_Control.'
179 else if (index(controlcard,'MCMA').gt.0) then
181 else if (index(controlcard,'MCEE').gt.0) then
183 else if (index(controlcard,'MULTCONF').gt.0) then
185 else if (index(controlcard,'MAP').gt.0) then
187 call readi(controlcard,'MAP',nmap,0)
188 else if (index(controlcard,'CSA').gt.0) then
190 crc else if (index(controlcard,'ZSCORE').gt.0) then
192 crc ZSCORE is rm from UNRES, modecalc=9 is available
195 cfcm else if (index(controlcard,'MCMF').gt.0) then
197 else if (index(controlcard,'SOFTREG').gt.0) then
199 else if (index(controlcard,'CHECK_BOND').gt.0) then
201 else if (index(controlcard,'TEST').gt.0) then
203 else if (index(controlcard,'MD').gt.0) then
205 else if (index(controlcard,'RE ').gt.0) then
209 lmuca=index(controlcard,'MUCA').gt.0
210 call readi(controlcard,'MUCADYN',mucadyn,0)
211 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
212 if (lmuca .and. (me.eq.king .or. .not.out1file ))
214 write (iout,*) 'MUCADYN=',mucadyn
215 write (iout,*) 'MUCASMOOTH=',muca_smooth
218 iscode=index(controlcard,'ONE_LETTER')
219 indphi=index(controlcard,'PHI')
220 indback=index(controlcard,'BACK')
221 iranconf=index(controlcard,'RAND_CONF')
222 i2ndstr=index(controlcard,'USE_SEC_PRED')
223 gradout=index(controlcard,'GRADOUT').gt.0
224 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
225 C DISTCHAINMAX become obsolete for periodic boundry condition
226 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
227 C Reading the dimensions of box in x,y,z coordinates
228 call reada(controlcard,'BOXX',boxxsize,100.0d0)
229 call reada(controlcard,'BOXY',boxysize,100.0d0)
230 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
231 c Cutoff range for interactions
232 call reada(controlcard,"R_CUT",r_cut,15.0d0)
233 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
234 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
235 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
236 if (lipthick.gt.0.0d0) then
237 bordliptop=(boxzsize+lipthick)/2.0
238 bordlipbot=bordliptop-lipthick
240 if ((bordliptop.gt.boxzsize).or.(borlipbot.lt.0.0))
241 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
242 buflipbot=bordlipbot+lipbufthick
243 bufliptop=bordliptop-lipbufthick
244 if ((lipbufthick*2.0d0).gt.lipthick)
245 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
247 write(iout,*) "bordliptop=",bordliptop
248 write(iout,*) "bordlipbot=",bordlipbot
249 write(iout,*) "bufliptop=",bufliptop
250 write(iout,*) "buflipbot=",buflipbot
253 if (me.eq.king .or. .not.out1file )
254 & write (iout,*) "DISTCHAINMAX",distchainmax
256 if(me.eq.king.or..not.out1file)
257 & write (iout,'(2a)') diagmeth(kdiag),
258 & ' routine used to diagonalize matrices.'
261 c--------------------------------------------------------------------------
262 subroutine read_REMDpar
266 implicit real*8 (a-h,o-z)
268 include 'COMMON.IOUNITS'
269 include 'COMMON.TIME1'
272 include 'COMMON.LANGEVIN'
274 include 'COMMON.LANGEVIN.lang0'
276 include 'COMMON.INTERACT'
277 include 'COMMON.NAMES'
279 include 'COMMON.REMD'
280 include 'COMMON.CONTROL'
281 include 'COMMON.SETUP'
283 character*320 controlcard
284 character*3200 controlcard1
285 integer iremd_m_total
287 if(me.eq.king.or..not.out1file)
288 & write (iout,*) "REMD setup"
290 call card_concat(controlcard)
291 call readi(controlcard,"NREP",nrep,3)
292 call readi(controlcard,"NSTEX",nstex,1000)
293 call reada(controlcard,"RETMIN",retmin,10.0d0)
294 call reada(controlcard,"RETMAX",retmax,1000.0d0)
295 mremdsync=(index(controlcard,'SYNC').gt.0)
296 call readi(controlcard,"NSYN",i_sync_step,100)
297 restart1file=(index(controlcard,'REST1FILE').gt.0)
298 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
299 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
300 if(max_cache_traj_use.gt.max_cache_traj)
301 & max_cache_traj_use=max_cache_traj
302 if(me.eq.king.or..not.out1file) then
303 cd if (traj1file) then
304 crc caching is in testing - NTWX is not ignored
305 cd write (iout,*) "NTWX value is ignored"
306 cd write (iout,*) " trajectory is stored to one file by master"
307 cd write (iout,*) " before exchange at NSTEX intervals"
309 write (iout,*) "NREP= ",nrep
310 write (iout,*) "NSTEX= ",nstex
311 write (iout,*) "SYNC= ",mremdsync
312 write (iout,*) "NSYN= ",i_sync_step
313 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
316 if (index(controlcard,'TLIST').gt.0) then
318 call card_concat(controlcard1)
319 read(controlcard1,*) (remd_t(i),i=1,nrep)
320 if(me.eq.king.or..not.out1file)
321 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
324 if (index(controlcard,'MLIST').gt.0) then
326 call card_concat(controlcard1)
327 read(controlcard1,*) (remd_m(i),i=1,nrep)
328 if(me.eq.king.or..not.out1file) then
329 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
332 iremd_m_total=iremd_m_total+remd_m(i)
334 write (iout,*) 'Total number of replicas ',iremd_m_total
337 if(me.eq.king.or..not.out1file)
338 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
341 c--------------------------------------------------------------------------
342 subroutine read_MDpar
346 implicit real*8 (a-h,o-z)
348 include 'COMMON.IOUNITS'
349 include 'COMMON.TIME1'
352 include 'COMMON.LANGEVIN'
354 include 'COMMON.LANGEVIN.lang0'
356 include 'COMMON.INTERACT'
357 include 'COMMON.NAMES'
359 include 'COMMON.SETUP'
360 include 'COMMON.CONTROL'
361 include 'COMMON.SPLITELE'
363 character*320 controlcard
365 call card_concat(controlcard)
366 call readi(controlcard,"NSTEP",n_timestep,1000000)
367 call readi(controlcard,"NTWE",ntwe,100)
368 call readi(controlcard,"NTWX",ntwx,1000)
369 call reada(controlcard,"DT",d_time,1.0d-1)
370 call reada(controlcard,"DVMAX",dvmax,2.0d1)
371 call reada(controlcard,"DAMAX",damax,1.0d1)
372 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
373 call readi(controlcard,"LANG",lang,0)
374 RESPA = index(controlcard,"RESPA") .gt. 0
375 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
376 ntime_split0=ntime_split
377 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
378 ntime_split0=ntime_split
379 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
380 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
381 rest = index(controlcard,"REST").gt.0
382 tbf = index(controlcard,"TBF").gt.0
383 usampl = index(controlcard,"USAMPL").gt.0
385 mdpdb = index(controlcard,"MDPDB").gt.0
386 call reada(controlcard,"T_BATH",t_bath,300.0d0)
387 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
388 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
389 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
390 if (count_reset_moment.eq.0) count_reset_moment=1000000000
391 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
392 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
393 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
394 if (count_reset_vel.eq.0) count_reset_vel=1000000000
395 large = index(controlcard,"LARGE").gt.0
396 print_compon = index(controlcard,"PRINT_COMPON").gt.0
397 rattle = index(controlcard,"RATTLE").gt.0
398 c if performing umbrella sampling, fragments constrained are read from the fragment file
404 if(me.eq.king.or..not.out1file) then
406 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
408 write (iout,'(a)') "The units are:"
409 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
410 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
411 & " acceleration: angstrom/(48.9 fs)**2"
412 write (iout,'(a)') "energy: kcal/mol, temperature: K"
414 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
415 write (iout,'(a60,f10.5,a)')
416 & "Initial time step of numerical integration:",d_time,
418 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
420 write (iout,'(2a,i4,a)')
421 & "A-MTS algorithm used; initial time step for fast-varying",
422 & " short-range forces split into",ntime_split," steps."
423 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
424 & r_cut," lambda",rlamb
426 write (iout,'(2a,f10.5)')
427 & "Maximum acceleration threshold to reduce the time step",
428 & "/increase split number:",damax
429 write (iout,'(2a,f10.5)')
430 & "Maximum predicted energy drift to reduce the timestep",
431 & "/increase split number:",edriftmax
432 write (iout,'(a60,f10.5)')
433 & "Maximum velocity threshold to reduce velocities:",dvmax
434 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
435 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
436 if (rattle) write (iout,'(a60)')
437 & "Rattle algorithm used to constrain the virtual bonds"
441 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
442 call reada(controlcard,"RWAT",rwat,1.4d0)
443 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
444 surfarea=index(controlcard,"SURFAREA").gt.0
445 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
446 if(me.eq.king.or..not.out1file)then
447 write (iout,'(/a,$)') "Langevin dynamics calculation"
450 & " with direct integration of Langevin equations"
451 else if (lang.eq.2) then
452 write (iout,'(a/)') " with TINKER stochasic MD integrator"
453 else if (lang.eq.3) then
454 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
455 else if (lang.eq.4) then
456 write (iout,'(a/)') " in overdamped mode"
458 write (iout,'(//a,i5)')
459 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
462 write (iout,'(a60,f10.5)') "Temperature:",t_bath
463 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
464 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
465 write (iout,'(a60,f10.5)')
466 & "Scaling factor of the friction forces:",scal_fric
467 if (surfarea) write (iout,'(2a,i10,a)')
468 & "Friction coefficients will be scaled by solvent-accessible",
469 & " surface area every",reset_fricmat," steps."
471 c Calculate friction coefficients and bounds of stochastic forces
472 eta=6*pi*cPoise*etawat
473 if(me.eq.king.or..not.out1file)
474 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
476 gamp=scal_fric*(pstok+rwat)*eta
477 stdfp=dsqrt(2*Rb*t_bath/d_time)
479 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
480 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
482 if(me.eq.king.or..not.out1file)then
483 write (iout,'(/2a/)')
484 & "Radii of site types and friction coefficients and std's of",
485 & " stochastic forces of fully exposed sites"
486 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
488 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
489 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
493 if(me.eq.king.or..not.out1file)then
494 write (iout,'(a)') "Berendsen bath calculation"
495 write (iout,'(a60,f10.5)') "Temperature:",t_bath
496 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
498 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
499 & count_reset_moment," steps"
501 & write (iout,'(a,i10,a)')
502 & "Velocities will be reset at random every",count_reset_vel,
506 if(me.eq.king.or..not.out1file)
507 & write (iout,'(a31)') "Microcanonical mode calculation"
509 if(me.eq.king.or..not.out1file)then
510 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
512 write(iout,*) "MD running with constraints."
513 write(iout,*) "Equilibration time ", eq_time, " mtus."
514 write(iout,*) "Constraining ", nfrag," fragments."
515 write(iout,*) "Length of each fragment, weight and q0:"
517 write (iout,*) "Set of restraints #",iset
519 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
520 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
522 write(iout,*) "constraints between ", npair, "fragments."
523 write(iout,*) "constraint pairs, weights and q0:"
525 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
526 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
528 write(iout,*) "angle constraints within ", nfrag_back,
529 & "backbone fragments."
530 write(iout,*) "fragment, weights:"
532 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
533 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
534 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
537 iset=mod(kolor,nset)+1
540 if(me.eq.king.or..not.out1file)
541 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
544 c------------------------------------------------------------------------------
547 C Read molecular data.
549 implicit real*8 (a-h,o-z)
555 include 'COMMON.IOUNITS'
558 include 'COMMON.INTERACT'
559 include 'COMMON.LOCAL'
560 include 'COMMON.NAMES'
561 include 'COMMON.CHAIN'
562 include 'COMMON.FFIELD'
563 include 'COMMON.SBRIDGE'
564 include 'COMMON.HEADER'
565 include 'COMMON.CONTROL'
566 include 'COMMON.DBASE'
567 include 'COMMON.THREAD'
568 include 'COMMON.CONTACTS'
569 include 'COMMON.TORCNSTR'
570 include 'COMMON.TIME1'
571 include 'COMMON.BOUNDS'
573 include 'COMMON.SETUP'
574 character*4 sequence(maxres)
576 double precision x(maxvar)
577 character*256 pdbfile
578 character*400 weightcard
579 character*80 weightcard_t,ucase
580 dimension itype_pdb(maxres)
581 common /pizda/ itype_pdb
582 logical seq_comp,fail
583 double precision energia(0:n_ene)
589 C Read weights of the subsequent energy terms.
590 call card_concat(weightcard)
591 call reada(weightcard,'WLONG',wlong,1.0D0)
592 call reada(weightcard,'WSC',wsc,wlong)
593 call reada(weightcard,'WSCP',wscp,wlong)
594 call reada(weightcard,'WELEC',welec,1.0D0)
595 call reada(weightcard,'WVDWPP',wvdwpp,welec)
596 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
597 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
598 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
599 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
600 call reada(weightcard,'WTURN3',wturn3,1.0D0)
601 call reada(weightcard,'WTURN4',wturn4,1.0D0)
602 call reada(weightcard,'WTURN6',wturn6,1.0D0)
603 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
604 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
605 call reada(weightcard,'WBOND',wbond,1.0D0)
606 call reada(weightcard,'WTOR',wtor,1.0D0)
607 call reada(weightcard,'WTORD',wtor_d,1.0D0)
608 call reada(weightcard,'WANG',wang,1.0D0)
609 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
610 call reada(weightcard,'SCAL14',scal14,0.4D0)
611 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
612 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
613 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
614 call reada(weightcard,'TEMP0',temp0,300.0d0)
615 call reada(weightcard,'WLT',wliptran,0.0D0)
616 if (index(weightcard,'SOFT').gt.0) ipot=6
617 C 12/1/95 Added weight for the multi-body term WCORR
618 call reada(weightcard,'WCORRH',wcorr,1.0D0)
619 if (wcorr4.gt.0.0d0) wcorr=wcorr4
639 if(me.eq.king.or..not.out1file)
640 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
641 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
643 10 format (/'Energy-term weights (unscaled):'//
644 & 'WSCC= ',f10.6,' (SC-SC)'/
645 & 'WSCP= ',f10.6,' (SC-p)'/
646 & 'WELEC= ',f10.6,' (p-p electr)'/
647 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
648 & 'WBOND= ',f10.6,' (stretching)'/
649 & 'WANG= ',f10.6,' (bending)'/
650 & 'WSCLOC= ',f10.6,' (SC local)'/
651 & 'WTOR= ',f10.6,' (torsional)'/
652 & 'WTORD= ',f10.6,' (double torsional)'/
653 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
654 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
655 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
656 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
657 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
658 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
659 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
660 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
661 & 'WTURN6= ',f10.6,' (turns, 6th order)')
662 if(me.eq.king.or..not.out1file)then
663 if (wcorr4.gt.0.0d0) then
664 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
665 & 'between contact pairs of peptide groups'
666 write (iout,'(2(a,f5.3/))')
667 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
668 & 'Range of quenching the correlation terms:',2*delt_corr
669 else if (wcorr.gt.0.0d0) then
670 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
671 & 'between contact pairs of peptide groups'
673 write (iout,'(a,f8.3)')
674 & 'Scaling factor of 1,4 SC-p interactions:',scal14
675 write (iout,'(a,f8.3)')
676 & 'General scaling factor of SC-p interactions:',scalscp
678 r0_corr=cutoff_corr-delt_corr
680 aad(i,1)=scalscp*aad(i,1)
681 aad(i,2)=scalscp*aad(i,2)
682 bad(i,1)=scalscp*bad(i,1)
683 bad(i,2)=scalscp*bad(i,2)
685 call rescale_weights(t_bath)
686 if(me.eq.king.or..not.out1file)
687 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
688 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
690 22 format (/'Energy-term weights (scaled):'//
691 & 'WSCC= ',f10.6,' (SC-SC)'/
692 & 'WSCP= ',f10.6,' (SC-p)'/
693 & 'WELEC= ',f10.6,' (p-p electr)'/
694 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
695 & 'WBOND= ',f10.6,' (stretching)'/
696 & 'WANG= ',f10.6,' (bending)'/
697 & 'WSCLOC= ',f10.6,' (SC local)'/
698 & 'WTOR= ',f10.6,' (torsional)'/
699 & 'WTORD= ',f10.6,' (double torsional)'/
700 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
701 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
702 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
703 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
704 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
705 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
706 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
707 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
708 & 'WTURN6= ',f10.6,' (turns, 6th order)')
709 if(me.eq.king.or..not.out1file)
710 & write (iout,*) "Reference temperature for weights calculation:",
712 call reada(weightcard,"D0CM",d0cm,3.78d0)
713 call reada(weightcard,"AKCM",akcm,15.1d0)
714 call reada(weightcard,"AKTH",akth,11.0d0)
715 call reada(weightcard,"AKCT",akct,12.0d0)
716 call reada(weightcard,"V1SS",v1ss,-1.08d0)
717 call reada(weightcard,"V2SS",v2ss,7.61d0)
718 call reada(weightcard,"V3SS",v3ss,13.7d0)
719 call reada(weightcard,"EBR",ebr,-5.50D0)
720 call reada(weightcard,"ATRISS",atriss,0.301D0)
721 call reada(weightcard,"BTRISS",btriss,0.021D0)
722 call reada(weightcard,"CTRISS",ctriss,1.001D0)
723 call reada(weightcard,"DTRISS",dtriss,1.001D0)
724 write (iout,*) "ATRISS=", atriss
725 write (iout,*) "BTRISS=", btriss
726 write (iout,*) "CTRISS=", ctriss
727 write (iout,*) "DTRISS=", dtriss
728 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
730 dyn_ss_mask(i)=.false.
734 dyn_ssbond_ij(i,j)=1.0d300
737 call reada(weightcard,"HT",Ht,0.0D0)
739 ss_depth=ebr/wsc-0.25*eps(1,1)
740 Ht=Ht/wsc-0.25*eps(1,1)
741 akcm=akcm*wstrain/wsc
742 akth=akth*wstrain/wsc
743 akct=akct*wstrain/wsc
744 v1ss=v1ss*wstrain/wsc
745 v2ss=v2ss*wstrain/wsc
746 v3ss=v3ss*wstrain/wsc
748 if (wstrain.ne.0.0) then
749 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
755 if(me.eq.king.or..not.out1file) then
756 write (iout,*) "Parameters of the SS-bond potential:"
757 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
759 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
760 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
761 write (iout,*)" HT",Ht
762 print *,'indpdb=',indpdb,' pdbref=',pdbref
764 if (indpdb.gt.0 .or. pdbref) then
765 read(inp,'(a)') pdbfile
766 if(me.eq.king.or..not.out1file)
767 & write (iout,'(2a)') 'PDB data will be read from file ',
768 & pdbfile(:ilen(pdbfile))
769 open(ipdbin,file=pdbfile,status='old',err=33)
771 33 write (iout,'(a)') 'Error opening PDB file.'
774 c write (iout,*) 'Begin reading pdb data'
777 c write (iout,*) 'Finished reading pdb data'
779 if(me.eq.king.or..not.out1file)
780 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
781 & ' nstart_sup=',nstart_sup
783 itype_pdb(i)=itype(i)
787 nct=nstart_sup+nsup-1
788 call contact(.false.,ncont_ref,icont_ref,co)
791 if(me.eq.king.or..not.out1file)
792 & write(iout,*)'Adding sidechains'
796 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
799 do while (fail.and.nsi.le.maxsi)
800 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
803 if(fail) write(iout,*)'Adding sidechain failed for res ',
804 & i,' after ',nsi,' trials'
809 if (indpdb.eq.0) then
810 C Read sequence if not taken from the pdb file.
812 c print *,'nres=',nres
813 if (iscode.gt.0) then
814 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
816 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
818 C Convert sequence to numeric code
820 itype(i)=rescode(i,sequence(i),iscode)
822 C Assign initial virtual bond lengths
828 vbld(i+nres)=dsc(iabs(itype(i)))
829 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
830 c write (iout,*) "i",i," itype",itype(i),
831 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
835 c print '(20i4)',(itype(i),i=1,nres)
838 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
840 if (itype(i).eq.ntyp1) then
844 else if (iabs(itype(i+1)).ne.20) then
846 else if (iabs(itype(i)).ne.20) then
853 if(me.eq.king.or..not.out1file)then
854 write (iout,*) "ITEL"
856 write (iout,*) i,itype(i),itel(i)
858 print *,'Call Read_Bridge.'
861 C 8/13/98 Set limits to generating the dihedral angles
866 read (inp,*) ndih_constr
867 if (ndih_constr.gt.0) then
869 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
871 if(me.eq.king.or..not.out1file)then
873 & 'There are',ndih_constr,' constraints on phi angles.'
875 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
880 phi0(i)=deg2rad*phi0(i)
881 drange(i)=deg2rad*drange(i)
883 C if(me.eq.king.or..not.out1file)
884 C & write (iout,*) 'FTORS',ftors
887 phibound(1,ii) = phi0(i)-drange(i)
888 phibound(2,ii) = phi0(i)+drange(i)
891 C first setting the theta boundaries to 0 to pi
892 C this mean that there is no energy penalty for any angle occuring this can be applied
893 C for generate random conformation but is not implemented in this way
898 C begin reading theta constrains this is quartic constrains allowing to
899 C have smooth second derivative
900 if (with_theta_constr) then
901 C with_theta_constr is keyword allowing for occurance of theta constrains
902 read (inp,*) ntheta_constr
903 C ntheta_constr is the number of theta constrains
904 if (ntheta_constr.gt.0) then
906 read (inp,*) (itheta_constr(i),theta_constr0(i),
907 & theta_drange(i),for_thet_constr(i),
909 C the above code reads from 1 to ntheta_constr
910 C itheta_constr(i) residue i for which is theta_constr
911 C theta_constr0 the global minimum value
912 C theta_drange is range for which there is no energy penalty
913 C for_thet_constr is the force constant for quartic energy penalty
915 if(me.eq.king.or..not.out1file)then
917 & 'There are',ntheta_constr,' constraints on phi angles.'
919 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
925 theta_constr0(i)=deg2rad*theta_constr0(i)
926 theta_drange(i)=deg2rad*theta_drange(i)
928 C if(me.eq.king.or..not.out1file)
929 C & write (iout,*) 'FTORS',ftors
930 C do i=1,ntheta_constr
931 C ii = itheta_constr(i)
932 C thetabound(1,ii) = phi0(i)-drange(i)
933 C thetabound(2,ii) = phi0(i)+drange(i)
935 endif ! ntheta_constr.gt.0
936 endif! with_theta_constr
938 C with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
939 C write (iout,*) "with_dihed_constr ",with_dihed_constr
944 write (iout,'(a)') 'Boundaries in phi angle sampling:'
946 write (iout,'(a3,i5,2f10.1)')
947 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
953 cd print *,'NNT=',NNT,' NCT=',NCT
954 if (itype(1).eq.ntyp1) nnt=2
955 if (itype(nres).eq.ntyp1) nct=nct-1
957 if(me.eq.king.or..not.out1file)
958 & write (iout,'(a,i3)') 'nsup=',nsup
960 if (nsup.le.(nct-nnt+1)) then
961 do i=0,nct-nnt+1-nsup
962 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
968 & 'Error - sequences to be superposed do not match.'
971 do i=0,nsup-(nct-nnt+1)
972 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
974 nstart_sup=nstart_sup+i
980 & 'Error - sequences to be superposed do not match.'
983 if (nsup.eq.0) nsup=nct-nnt
984 if (nstart_sup.eq.0) nstart_sup=nnt
985 if (nstart_seq.eq.0) nstart_seq=nnt
986 if(me.eq.king.or..not.out1file)
987 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
988 & ' nstart_seq=',nstart_seq
990 c--- Zscore rms -------
991 if (nz_start.eq.0) nz_start=nnt
992 if (nz_end.eq.0 .and. nsup.gt.0) then
994 else if (nz_end.eq.0) then
997 if(me.eq.king.or..not.out1file)then
998 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
999 write (iout,*) 'IZ_SC=',iz_sc
1001 c----------------------
1004 if (.not.pdbref) then
1005 call read_angles(inp,*38)
1007 38 write (iout,'(a)') 'Error reading reference structure.'
1009 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1010 stop 'Error reading reference structure'
1014 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1021 cref(j,i,kkk)=c(j,i)
1024 call contact(.true.,ncont_ref,icont_ref,co)
1028 write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1030 if (constr_dist.gt.0) call read_dist_constr
1031 write (iout,*) "After read_dist_constr nhpb",nhpb
1032 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1034 if(me.eq.king.or..not.out1file)
1035 & write (iout,*) 'Contact order:',co
1037 if(me.eq.king.or..not.out1file)
1038 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1041 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1043 if(me.eq.king.or..not.out1file)
1044 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1045 & icont_ref(1,i),' ',
1046 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1050 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1051 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1052 & modecalc.ne.10) then
1053 C If input structure hasn't been supplied from the PDB file read or generate
1055 if (iranconf.eq.0 .and. .not. extconf) then
1056 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1057 & write (iout,'(a)') 'Initial geometry will be read in.'
1059 read(inp,'(8f10.5)',end=36,err=36)
1060 & ((c(l,k),l=1,3),k=1,nres),
1061 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1062 write (iout,*) "Exit READ_CART"
1063 write (iout,'(8f10.5)')
1064 & ((c(l,k),l=1,3),k=1,nres),
1065 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1066 call int_from_cart1(.true.)
1067 write (iout,*) "Finish INT_TO_CART"
1070 dc(j,i)=c(j,i+1)-c(j,i)
1071 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1075 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1077 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1078 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1084 call read_angles(inp,*36)
1087 36 write (iout,'(a)') 'Error reading angle file.'
1089 call mpi_finalize( MPI_COMM_WORLD,IERR )
1091 stop 'Error reading angle file.'
1093 else if (extconf) then
1094 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1095 & write (iout,'(a)') 'Extended chain initial geometry.'
1097 theta(i)=90d0*deg2rad
1100 phi(i)=180d0*deg2rad
1103 alph(i)=110d0*deg2rad
1106 omeg(i)=-120d0*deg2rad
1107 if (itype(i).le.0) omeg(i)=-omeg(i)
1110 if(me.eq.king.or..not.out1file)
1111 & write (iout,'(a)') 'Random-generated initial geometry.'
1115 if (me.eq.king .or. fg_rank.eq.0 .and. (
1116 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1120 call gen_rand_conf(itmp,*30)
1122 30 write (iout,*) 'Failed to generate random conformation',
1123 & ', itrial=',itrial
1124 write (*,*) 'Processor:',me,
1125 & ' Failed to generate random conformation',
1135 write (iout,'(a,i3,a)') 'Processor:',me,
1136 & ' error in generating random conformation.'
1137 write (*,'(a,i3,a)') 'Processor:',me,
1138 & ' error in generating random conformation.'
1141 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1146 & ' error in generating random conformation.'
1151 elseif (modecalc.eq.4) then
1152 read (inp,'(a)') intinname
1153 open (intin,file=intinname,status='old',err=333)
1154 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1155 & write (iout,'(a)') 'intinname',intinname
1156 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1158 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1160 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1162 stop 'Error opening angle file.'
1166 C Generate distance constraints, if the PDB structure is to be regularized.
1167 if (nthread.gt.0) then
1168 call read_threadbase
1171 if (me.eq.king .or. .not. out1file)
1173 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1174 write (iout,'(/a,i3,a)')
1175 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1176 write (iout,'(20i4)') (iss(i),i=1,ns)
1178 write(iout,*)"Running with dynamic disulfide-bond formation"
1180 write (iout,'(/a/)') 'Pre-formed links are:'
1186 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1187 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1193 if (ns.gt.0.and.dyn_ss) then
1197 forcon(i-nss)=forcon(i)
1204 dyn_ss_mask(iss(i))=.true.
1207 if (i2ndstr.gt.0) call secstrp2dihc
1208 c call geom_to_var(nvar,x)
1209 c call etotal(energia(0))
1210 c call enerprint(energia(0))
1211 c call briefout(0,etot)
1213 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1214 cd write (iout,'(a)') 'Variable list:'
1215 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1217 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1218 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1219 & 'Processor',myrank,': end reading molecular data.'
1224 c--------------------------------------------------------------------------
1225 logical function seq_comp(itypea,itypeb,length)
1227 integer length,itypea(length),itypeb(length)
1230 if (itypea(i).ne.itypeb(i)) then
1238 c-----------------------------------------------------------------------------
1239 subroutine read_bridge
1240 C Read information about disulfide bridges.
1241 implicit real*8 (a-h,o-z)
1242 include 'DIMENSIONS'
1246 include 'COMMON.IOUNITS'
1247 include 'COMMON.GEO'
1248 include 'COMMON.VAR'
1249 include 'COMMON.INTERACT'
1250 include 'COMMON.LOCAL'
1251 include 'COMMON.NAMES'
1252 include 'COMMON.CHAIN'
1253 include 'COMMON.FFIELD'
1254 include 'COMMON.SBRIDGE'
1255 include 'COMMON.HEADER'
1256 include 'COMMON.CONTROL'
1257 include 'COMMON.DBASE'
1258 include 'COMMON.THREAD'
1259 include 'COMMON.TIME1'
1260 include 'COMMON.SETUP'
1261 C Read bridging residues.
1262 read (inp,*) ns,(iss(i),i=1,ns)
1264 if(me.eq.king.or..not.out1file)
1265 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1266 C Check whether the specified bridging residues are cystines.
1268 if (itype(iss(i)).ne.1) then
1269 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1270 & 'Do you REALLY think that the residue ',
1271 & restyp(itype(iss(i))),i,
1272 & ' can form a disulfide bridge?!!!'
1273 write (*,'(2a,i3,a)')
1274 & 'Do you REALLY think that the residue ',
1275 & restyp(itype(iss(i))),i,
1276 & ' can form a disulfide bridge?!!!'
1278 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1283 C Read preformed bridges.
1285 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1287 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1290 C Check if the residues involved in bridges are in the specified list of
1291 C bridging residues.
1294 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1295 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1296 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1297 & ' contains residues present in other pairs.'
1298 write (*,'(a,i3,a)') 'Disulfide pair',i,
1299 & ' contains residues present in other pairs.'
1301 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1307 if (ihpb(i).eq.iss(j)) goto 10
1309 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1312 if (jhpb(i).eq.iss(j)) goto 20
1314 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1320 ihpb(i)=ihpb(i)+nres
1321 jhpb(i)=jhpb(i)+nres
1327 c----------------------------------------------------------------------------
1328 subroutine read_x(kanal,*)
1329 implicit real*8 (a-h,o-z)
1330 include 'DIMENSIONS'
1331 include 'COMMON.GEO'
1332 include 'COMMON.VAR'
1333 include 'COMMON.CHAIN'
1334 include 'COMMON.IOUNITS'
1335 include 'COMMON.CONTROL'
1336 include 'COMMON.LOCAL'
1337 include 'COMMON.INTERACT'
1338 c Read coordinates from input
1340 read(kanal,'(8f10.5)',end=10,err=10)
1341 & ((c(l,k),l=1,3),k=1,nres),
1342 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1345 c(j,2*nres)=c(j,nres)
1347 call int_from_cart1(.false.)
1350 dc(j,i)=c(j,i+1)-c(j,i)
1351 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1355 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1357 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1358 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1366 c----------------------------------------------------------------------------
1367 subroutine read_threadbase
1368 implicit real*8 (a-h,o-z)
1369 include 'DIMENSIONS'
1370 include 'COMMON.IOUNITS'
1371 include 'COMMON.GEO'
1372 include 'COMMON.VAR'
1373 include 'COMMON.INTERACT'
1374 include 'COMMON.LOCAL'
1375 include 'COMMON.NAMES'
1376 include 'COMMON.CHAIN'
1377 include 'COMMON.FFIELD'
1378 include 'COMMON.SBRIDGE'
1379 include 'COMMON.HEADER'
1380 include 'COMMON.CONTROL'
1381 include 'COMMON.DBASE'
1382 include 'COMMON.THREAD'
1383 include 'COMMON.TIME1'
1384 C Read pattern database for threading.
1385 read (icbase,*) nseq
1387 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1388 & nres_base(2,i),nres_base(3,i)
1389 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1391 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1392 c & nres_base(2,i),nres_base(3,i)
1393 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1397 if (weidis.eq.0.0D0) weidis=0.1D0
1406 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1407 write (iout,'(a,i5)') 'nexcl: ',nexcl
1408 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1411 c------------------------------------------------------------------------------
1412 subroutine setup_var
1413 implicit real*8 (a-h,o-z)
1414 include 'DIMENSIONS'
1415 include 'COMMON.IOUNITS'
1416 include 'COMMON.GEO'
1417 include 'COMMON.VAR'
1418 include 'COMMON.INTERACT'
1419 include 'COMMON.LOCAL'
1420 include 'COMMON.NAMES'
1421 include 'COMMON.CHAIN'
1422 include 'COMMON.FFIELD'
1423 include 'COMMON.SBRIDGE'
1424 include 'COMMON.HEADER'
1425 include 'COMMON.CONTROL'
1426 include 'COMMON.DBASE'
1427 include 'COMMON.THREAD'
1428 include 'COMMON.TIME1'
1429 C Set up variable list.
1435 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1437 ialph(i,1)=nvar+nside
1441 if (indphi.gt.0) then
1443 else if (indback.gt.0) then
1448 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1451 c----------------------------------------------------------------------------
1452 subroutine gen_dist_constr
1453 C Generate CA distance constraints.
1454 implicit real*8 (a-h,o-z)
1455 include 'DIMENSIONS'
1456 include 'COMMON.IOUNITS'
1457 include 'COMMON.GEO'
1458 include 'COMMON.VAR'
1459 include 'COMMON.INTERACT'
1460 include 'COMMON.LOCAL'
1461 include 'COMMON.NAMES'
1462 include 'COMMON.CHAIN'
1463 include 'COMMON.FFIELD'
1464 include 'COMMON.SBRIDGE'
1465 include 'COMMON.HEADER'
1466 include 'COMMON.CONTROL'
1467 include 'COMMON.DBASE'
1468 include 'COMMON.THREAD'
1469 include 'COMMON.TIME1'
1470 dimension itype_pdb(maxres)
1471 common /pizda/ itype_pdb
1473 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1474 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1475 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1477 do i=nstart_sup,nstart_sup+nsup-1
1478 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1479 cd & ' seq_pdb', restyp(itype_pdb(i))
1480 do j=i+2,nstart_sup+nsup-1
1482 ihpb(nhpb)=i+nstart_seq-nstart_sup
1483 jhpb(nhpb)=j+nstart_seq-nstart_sup
1485 dhpb(nhpb)=dist(i,j)
1488 cd write (iout,'(a)') 'Distance constraints:'
1493 cd if (ii.gt.nres) then
1498 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1499 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1500 cd & dhpb(i),forcon(i)
1504 c----------------------------------------------------------------------------
1506 implicit real*8 (a-h,o-z)
1507 include 'DIMENSIONS'
1508 include 'COMMON.MAP'
1509 include 'COMMON.IOUNITS'
1510 character*3 angid(4) /'THE','PHI','ALP','OME'/
1511 character*80 mapcard,ucase
1513 read (inp,'(a)') mapcard
1514 mapcard=ucase(mapcard)
1515 if (index(mapcard,'PHI').gt.0) then
1517 else if (index(mapcard,'THE').gt.0) then
1519 else if (index(mapcard,'ALP').gt.0) then
1521 else if (index(mapcard,'OME').gt.0) then
1524 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1525 stop 'Error - illegal variable spec in MAP card.'
1527 call readi (mapcard,'RES1',res1(imap),0)
1528 call readi (mapcard,'RES2',res2(imap),0)
1529 if (res1(imap).eq.0) then
1530 res1(imap)=res2(imap)
1531 else if (res2(imap).eq.0) then
1532 res2(imap)=res1(imap)
1534 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1536 & 'Error - illegal definition of variable group in MAP.'
1537 stop 'Error - illegal definition of variable group in MAP.'
1539 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1540 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1541 call readi(mapcard,'NSTEP',nstep(imap),0)
1542 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1544 & 'Illegal boundary and/or step size specification in MAP.'
1545 stop 'Illegal boundary and/or step size specification in MAP.'
1550 c----------------------------------------------------------------------------
1552 implicit real*8 (a-h,o-z)
1553 include 'DIMENSIONS'
1554 include 'COMMON.IOUNITS'
1555 include 'COMMON.GEO'
1556 include 'COMMON.CSA'
1557 include 'COMMON.BANK'
1558 include 'COMMON.CONTROL'
1560 character*620 mcmcard
1561 call card_concat(mcmcard)
1563 call readi(mcmcard,'NCONF',nconf,50)
1564 call readi(mcmcard,'NADD',nadd,0)
1565 call readi(mcmcard,'JSTART',jstart,1)
1566 call readi(mcmcard,'JEND',jend,1)
1567 call readi(mcmcard,'NSTMAX',nstmax,500000)
1568 call readi(mcmcard,'N0',n0,1)
1569 call readi(mcmcard,'N1',n1,6)
1570 call readi(mcmcard,'N2',n2,4)
1571 call readi(mcmcard,'N3',n3,0)
1572 call readi(mcmcard,'N4',n4,0)
1573 call readi(mcmcard,'N5',n5,0)
1574 call readi(mcmcard,'N6',n6,10)
1575 call readi(mcmcard,'N7',n7,0)
1576 call readi(mcmcard,'N8',n8,0)
1577 call readi(mcmcard,'N9',n9,0)
1578 call readi(mcmcard,'N14',n14,0)
1579 call readi(mcmcard,'N15',n15,0)
1580 call readi(mcmcard,'N16',n16,0)
1581 call readi(mcmcard,'N17',n17,0)
1582 call readi(mcmcard,'N18',n18,0)
1584 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1586 call readi(mcmcard,'NDIFF',ndiff,2)
1587 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1588 call readi(mcmcard,'IS1',is1,1)
1589 call readi(mcmcard,'IS2',is2,8)
1590 call readi(mcmcard,'NRAN0',nran0,4)
1591 call readi(mcmcard,'NRAN1',nran1,2)
1592 call readi(mcmcard,'IRR',irr,1)
1593 call readi(mcmcard,'NSEED',nseed,20)
1594 call readi(mcmcard,'NTOTAL',ntotal,10000)
1595 call reada(mcmcard,'CUT1',cut1,2.0d0)
1596 call reada(mcmcard,'CUT2',cut2,5.0d0)
1597 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1598 call readi(mcmcard,'ICMAX',icmax,3)
1599 call readi(mcmcard,'IRESTART',irestart,0)
1600 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1603 call reada(mcmcard,'DELE',dele,20.0d0)
1604 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1605 call readi(mcmcard,'IREF',iref,0)
1606 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1607 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1608 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1609 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1610 write (iout,*) "NCONF_IN",nconf_in
1613 c----------------------------------------------------------------------------
1614 cfmc subroutine mcmfread
1615 cfmc implicit real*8 (a-h,o-z)
1616 cfmc include 'DIMENSIONS'
1617 cfmc include 'COMMON.MCMF'
1618 cfmc include 'COMMON.IOUNITS'
1619 cfmc include 'COMMON.GEO'
1620 cfmc character*80 ucase
1621 cfmc character*620 mcmcard
1622 cfmc call card_concat(mcmcard)
1624 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1625 cfmc write(iout,*)'MAXRANT=',maxrant
1626 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1627 cfmc write(iout,*)'MAXFAM=',maxfam
1628 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1629 cfmc write(iout,*)'NNET1=',nnet1
1630 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1631 cfmc write(iout,*)'NNET2=',nnet2
1632 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1633 cfmc write(iout,*)'NNET3=',nnet3
1634 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1635 cfmc write(iout,*)'ILASTT=',ilastt
1636 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1637 cfmc write(iout,*)'MAXSTR=',maxstr
1638 cfmc maxstr_f=maxstr/maxfam
1639 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1640 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1641 cfmc write(iout,*)'NMCMF=',nmcmf
1642 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1643 cfmc write(iout,*)'IFOCUS=',ifocus
1644 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1645 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1646 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1647 cfmc write(iout,*)'INTPRT=',intprt
1648 cfmc call readi(mcmcard,'IPRT',iprt,100)
1649 cfmc write(iout,*)'IPRT=',iprt
1650 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1651 cfmc write(iout,*)'IMAXTR=',imaxtr
1652 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1653 cfmc write(iout,*)'MAXEVEN=',maxeven
1654 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1655 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1656 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1657 cfmc write(iout,*)'INIMIN=',inimin
1658 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1659 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1660 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1661 cfmc write(iout,*)'NTHREAD=',nthread
1662 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1663 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1664 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1665 cfmc write(iout,*)'MAXPERT=',maxpert
1666 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1667 cfmc write(iout,*)'IRMSD=',irmsd
1668 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1669 cfmc write(iout,*)'DENEMIN=',denemin
1670 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1671 cfmc write(iout,*)'RCUT1S=',rcut1s
1672 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1673 cfmc write(iout,*)'RCUT1E=',rcut1e
1674 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1675 cfmc write(iout,*)'RCUT2S=',rcut2s
1676 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1677 cfmc write(iout,*)'RCUT2E=',rcut2e
1678 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1679 cfmc write(iout,*)'DPERT1=',d_pert1
1680 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1681 cfmc write(iout,*)'DPERT1A=',d_pert1a
1682 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1683 cfmc write(iout,*)'DPERT2=',d_pert2
1684 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1685 cfmc write(iout,*)'DPERT2A=',d_pert2a
1686 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1687 cfmc write(iout,*)'DPERT2B=',d_pert2b
1688 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1689 cfmc write(iout,*)'DPERT2C=',d_pert2c
1690 cfmc d_pert1=deg2rad*d_pert1
1691 cfmc d_pert1a=deg2rad*d_pert1a
1692 cfmc d_pert2=deg2rad*d_pert2
1693 cfmc d_pert2a=deg2rad*d_pert2a
1694 cfmc d_pert2b=deg2rad*d_pert2b
1695 cfmc d_pert2c=deg2rad*d_pert2c
1696 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1697 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1698 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1699 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1700 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1701 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1702 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1703 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1704 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1705 cfmc write(iout,*)'RCUTINI=',rcutini
1706 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1707 cfmc write(iout,*)'GRAT=',grat
1708 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1709 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1713 c----------------------------------------------------------------------------
1715 implicit real*8 (a-h,o-z)
1716 include 'DIMENSIONS'
1717 include 'COMMON.MCM'
1718 include 'COMMON.MCE'
1719 include 'COMMON.IOUNITS'
1721 character*320 mcmcard
1722 call card_concat(mcmcard)
1723 call readi(mcmcard,'MAXACC',maxacc,100)
1724 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1725 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1726 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1727 call readi(mcmcard,'MAXREPM',maxrepm,200)
1728 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1729 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1730 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1731 call reada(mcmcard,'E_UP',e_up,5.0D0)
1732 call reada(mcmcard,'DELTE',delte,0.1D0)
1733 call readi(mcmcard,'NSWEEP',nsweep,5)
1734 call readi(mcmcard,'NSTEPH',nsteph,0)
1735 call readi(mcmcard,'NSTEPC',nstepc,0)
1736 call reada(mcmcard,'TMIN',tmin,298.0D0)
1737 call reada(mcmcard,'TMAX',tmax,298.0D0)
1738 call readi(mcmcard,'NWINDOW',nwindow,0)
1739 call readi(mcmcard,'PRINT_MC',print_mc,0)
1740 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1741 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1742 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1743 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1744 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1745 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1746 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1747 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1748 if (nwindow.gt.0) then
1749 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1751 winlen(i)=winend(i)-winstart(i)+1
1754 if (tmax.lt.tmin) tmax=tmin
1755 if (tmax.eq.tmin) then
1759 if (nstepc.gt.0 .and. nsteph.gt.0) then
1760 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1761 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1763 C Probabilities of different move types
1764 sumpro_type(0)=0.0D0
1765 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1766 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1767 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1768 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1769 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1770 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1771 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1773 print *,'i',i,' sumprotype',sumpro_type(i)
1774 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1775 print *,'i',i,' sumprotype',sumpro_type(i)
1779 c----------------------------------------------------------------------------
1780 subroutine read_minim
1781 implicit real*8 (a-h,o-z)
1782 include 'DIMENSIONS'
1783 include 'COMMON.MINIM'
1784 include 'COMMON.IOUNITS'
1786 character*320 minimcard
1787 call card_concat(minimcard)
1788 call readi(minimcard,'MAXMIN',maxmin,2000)
1789 call readi(minimcard,'MAXFUN',maxfun,5000)
1790 call readi(minimcard,'MINMIN',minmin,maxmin)
1791 call readi(minimcard,'MINFUN',minfun,maxmin)
1792 call reada(minimcard,'TOLF',tolf,1.0D-2)
1793 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1794 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1795 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1796 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1797 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1798 & 'Options in energy minimization:'
1799 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1800 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1801 & 'MinMin:',MinMin,' MinFun:',MinFun,
1802 & ' TolF:',TolF,' RTolF:',RTolF
1805 c----------------------------------------------------------------------------
1806 subroutine read_angles(kanal,*)
1807 implicit real*8 (a-h,o-z)
1808 include 'DIMENSIONS'
1809 include 'COMMON.GEO'
1810 include 'COMMON.VAR'
1811 include 'COMMON.CHAIN'
1812 include 'COMMON.IOUNITS'
1813 include 'COMMON.CONTROL'
1814 c Read angles from input
1816 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1817 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1818 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1819 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1822 c 9/7/01 avoid 180 deg valence angle
1823 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1825 theta(i)=deg2rad*theta(i)
1826 phi(i)=deg2rad*phi(i)
1827 alph(i)=deg2rad*alph(i)
1828 omeg(i)=deg2rad*omeg(i)
1833 c----------------------------------------------------------------------------
1834 subroutine reada(rekord,lancuch,wartosc,default)
1836 character*(*) rekord,lancuch
1837 double precision wartosc,default
1840 iread=index(rekord,lancuch)
1841 if (iread.eq.0) then
1845 iread=iread+ilen(lancuch)+1
1846 read (rekord(iread:),*,err=10,end=10) wartosc
1851 c----------------------------------------------------------------------------
1852 subroutine readi(rekord,lancuch,wartosc,default)
1854 character*(*) rekord,lancuch
1855 integer wartosc,default
1858 iread=index(rekord,lancuch)
1859 if (iread.eq.0) then
1863 iread=iread+ilen(lancuch)+1
1864 read (rekord(iread:),*,err=10,end=10) wartosc
1869 c----------------------------------------------------------------------------
1870 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1873 integer tablica(dim),default
1874 character*(*) rekord,lancuch
1881 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1882 if (iread.eq.0) return
1883 iread=iread+ilen(lancuch)+1
1884 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1887 c----------------------------------------------------------------------------
1888 subroutine multreada(rekord,lancuch,tablica,dim,default)
1891 double precision tablica(dim),default
1892 character*(*) rekord,lancuch
1899 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1900 if (iread.eq.0) return
1901 iread=iread+ilen(lancuch)+1
1902 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1905 c----------------------------------------------------------------------------
1906 subroutine openunits
1907 implicit real*8 (a-h,o-z)
1908 include 'DIMENSIONS'
1911 character*16 form,nodename
1914 include 'COMMON.SETUP'
1915 include 'COMMON.IOUNITS'
1917 include 'COMMON.CONTROL'
1918 integer lenpre,lenpot,ilen,lentmp
1920 character*3 out1file_text,ucase
1923 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1924 call getenv_loc("PREFIX",prefix)
1926 call getenv_loc("POT",pot)
1927 call getenv_loc("DIRTMP",tmpdir)
1928 call getenv_loc("CURDIR",curdir)
1929 call getenv_loc("OUT1FILE",out1file_text)
1930 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1931 out1file_text=ucase(out1file_text)
1932 if (out1file_text(1:1).eq."Y") then
1935 out1file=fg_rank.gt.0
1940 if (lentmp.gt.0) then
1941 write (*,'(80(1h!))')
1942 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1943 write (*,'(80(1h!))')
1944 write (*,*)"All output files will be on node /tmp directory."
1946 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1947 if (me.eq.king) then
1948 write (*,*) "The master node is ",nodename
1949 else if (fg_rank.eq.0) then
1950 write (*,*) "I am the CG slave node ",nodename
1952 write (*,*) "I am the FG slave node ",nodename
1955 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1956 lenpre = lentmp+lenpre+1
1958 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1959 C Get the names and open the input files
1960 #if defined(WINIFL) || defined(WINPGI)
1961 open(1,file=pref_orig(:ilen(pref_orig))//
1962 & '.inp',status='old',readonly,shared)
1963 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1964 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1965 C Get parameter filenames and open the parameter files.
1966 call getenv_loc('BONDPAR',bondname)
1967 open (ibond,file=bondname,status='old',readonly,shared)
1968 call getenv_loc('THETPAR',thetname)
1969 open (ithep,file=thetname,status='old',readonly,shared)
1970 call getenv_loc('ROTPAR',rotname)
1971 open (irotam,file=rotname,status='old',readonly,shared)
1972 call getenv_loc('TORPAR',torname)
1973 open (itorp,file=torname,status='old',readonly,shared)
1974 call getenv_loc('TORDPAR',tordname)
1975 open (itordp,file=tordname,status='old',readonly,shared)
1976 call getenv_loc('FOURIER',fouriername)
1977 open (ifourier,file=fouriername,status='old',readonly,shared)
1978 call getenv_loc('ELEPAR',elename)
1979 open (ielep,file=elename,status='old',readonly,shared)
1980 call getenv_loc('SIDEPAR',sidename)
1981 open (isidep,file=sidename,status='old',readonly,shared)
1982 #elif (defined CRAY) || (defined AIX)
1983 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1985 c print *,"Processor",myrank," opened file 1"
1986 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1987 c print *,"Processor",myrank," opened file 9"
1988 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1989 C Get parameter filenames and open the parameter files.
1990 call getenv_loc('BONDPAR',bondname)
1991 open (ibond,file=bondname,status='old',action='read')
1992 c print *,"Processor",myrank," opened file IBOND"
1993 call getenv_loc('THETPAR',thetname)
1994 open (ithep,file=thetname,status='old',action='read')
1995 c print *,"Processor",myrank," opened file ITHEP"
1996 call getenv_loc('ROTPAR',rotname)
1997 open (irotam,file=rotname,status='old',action='read')
1998 c print *,"Processor",myrank," opened file IROTAM"
1999 call getenv_loc('TORPAR',torname)
2000 open (itorp,file=torname,status='old',action='read')
2001 c print *,"Processor",myrank," opened file ITORP"
2002 call getenv_loc('TORDPAR',tordname)
2003 open (itordp,file=tordname,status='old',action='read')
2004 c print *,"Processor",myrank," opened file ITORDP"
2005 call getenv_loc('SCCORPAR',sccorname)
2006 open (isccor,file=sccorname,status='old',action='read')
2007 c print *,"Processor",myrank," opened file ISCCOR"
2008 call getenv_loc('FOURIER',fouriername)
2009 open (ifourier,file=fouriername,status='old',action='read')
2010 c print *,"Processor",myrank," opened file IFOURIER"
2011 call getenv_loc('ELEPAR',elename)
2012 open (ielep,file=elename,status='old',action='read')
2013 c print *,"Processor",myrank," opened file IELEP"
2014 call getenv_loc('SIDEPAR',sidename)
2015 open (isidep,file=sidename,status='old',action='read')
2016 c print *,"Processor",myrank," opened file ISIDEP"
2017 c print *,"Processor",myrank," opened parameter files"
2019 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2020 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2021 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2022 C Get parameter filenames and open the parameter files.
2023 call getenv_loc('BONDPAR',bondname)
2024 open (ibond,file=bondname,status='old')
2025 call getenv_loc('THETPAR',thetname)
2026 open (ithep,file=thetname,status='old')
2027 call getenv_loc('ROTPAR',rotname)
2028 open (irotam,file=rotname,status='old')
2029 call getenv_loc('TORPAR',torname)
2030 open (itorp,file=torname,status='old')
2031 call getenv_loc('TORDPAR',tordname)
2032 open (itordp,file=tordname,status='old')
2033 call getenv_loc('SCCORPAR',sccorname)
2034 open (isccor,file=sccorname,status='old')
2035 call getenv_loc('FOURIER',fouriername)
2036 open (ifourier,file=fouriername,status='old')
2037 call getenv_loc('ELEPAR',elename)
2038 open (ielep,file=elename,status='old')
2039 call getenv_loc('SIDEPAR',sidename)
2040 open (isidep,file=sidename,status='old')
2042 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2044 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2045 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2046 C Get parameter filenames and open the parameter files.
2047 call getenv_loc('BONDPAR',bondname)
2048 open (ibond,file=bondname,status='old',readonly)
2049 call getenv_loc('THETPAR',thetname)
2050 open (ithep,file=thetname,status='old',readonly)
2051 call getenv_loc('ROTPAR',rotname)
2052 open (irotam,file=rotname,status='old',readonly)
2053 call getenv_loc('TORPAR',torname)
2054 open (itorp,file=torname,status='old',readonly)
2055 call getenv_loc('TORDPAR',tordname)
2056 open (itordp,file=tordname,status='old',readonly)
2057 call getenv_loc('SCCORPAR',sccorname)
2058 open (isccor,file=sccorname,status='old',readonly)
2060 call getenv_loc('THETPARPDB',thetname_pdb)
2061 print *,"thetname_pdb ",thetname_pdb
2062 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2063 print *,ithep_pdb," opened"
2065 call getenv_loc('FOURIER',fouriername)
2066 open (ifourier,file=fouriername,status='old',readonly)
2067 call getenv_loc('ELEPAR',elename)
2068 open (ielep,file=elename,status='old',readonly)
2069 call getenv_loc('SIDEPAR',sidename)
2070 open (isidep,file=sidename,status='old',readonly)
2071 call getenv_loc('LIPTRANPAR',liptranname)
2072 open (iliptranpar,file=liptranname,status='old',action='read')
2074 call getenv_loc('ROTPARPDB',rotname_pdb)
2075 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2080 C 8/9/01 In the newest version SCp interaction constants are read from a file
2081 C Use -DOLDSCP to use hard-coded constants instead.
2083 call getenv_loc('SCPPAR',scpname)
2084 #if defined(WINIFL) || defined(WINPGI)
2085 open (iscpp,file=scpname,status='old',readonly,shared)
2086 #elif (defined CRAY) || (defined AIX)
2087 open (iscpp,file=scpname,status='old',action='read')
2089 open (iscpp,file=scpname,status='old')
2091 open (iscpp,file=scpname,status='old',readonly)
2094 call getenv_loc('PATTERN',patname)
2095 #if defined(WINIFL) || defined(WINPGI)
2096 open (icbase,file=patname,status='old',readonly,shared)
2097 #elif (defined CRAY) || (defined AIX)
2098 open (icbase,file=patname,status='old',action='read')
2100 open (icbase,file=patname,status='old')
2102 open (icbase,file=patname,status='old',readonly)
2105 C Open output file only for CG processes
2106 c print *,"Processor",myrank," fg_rank",fg_rank
2107 if (fg_rank.eq.0) then
2109 if (nodes.eq.1) then
2112 npos = dlog10(dfloat(nodes-1))+1
2114 if (npos.lt.3) npos=3
2115 write (liczba,'(i1)') npos
2116 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2118 write (liczba,form) me
2119 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2120 & liczba(:ilen(liczba))
2121 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2123 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2125 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2126 & liczba(:ilen(liczba))//'.mol2'
2127 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2128 & liczba(:ilen(liczba))//'.stat'
2130 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2131 & //liczba(:ilen(liczba))//'.stat')
2132 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2135 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2136 & liczba(:ilen(liczba))//'.const'
2141 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2142 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2143 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2144 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2145 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2147 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2149 rest2name=prefix(:ilen(prefix))//'.rst'
2151 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2154 #if defined(AIX) || defined(PGI)
2155 if (me.eq.king .or. .not. out1file)
2156 & open(iout,file=outname,status='unknown')
2158 if (fg_rank.gt.0) then
2159 write (liczba,'(i3.3)') myrank/nfgtasks
2160 write (ll,'(bz,i3.3)') fg_rank
2161 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2166 open(igeom,file=intname,status='unknown',position='append')
2167 open(ipdb,file=pdbname,status='unknown')
2168 open(imol2,file=mol2name,status='unknown')
2169 open(istat,file=statname,status='unknown',position='append')
2171 c1out open(iout,file=outname,status='unknown')
2174 if (me.eq.king .or. .not.out1file)
2175 & open(iout,file=outname,status='unknown')
2177 if (fg_rank.gt.0) then
2178 write (liczba,'(i3.3)') myrank/nfgtasks
2179 write (ll,'(bz,i3.3)') fg_rank
2180 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2185 open(igeom,file=intname,status='unknown',access='append')
2186 open(ipdb,file=pdbname,status='unknown')
2187 open(imol2,file=mol2name,status='unknown')
2188 open(istat,file=statname,status='unknown',access='append')
2190 c1out open(iout,file=outname,status='unknown')
2193 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2194 csa_seed=prefix(:lenpre)//'.CSA.seed'
2195 csa_history=prefix(:lenpre)//'.CSA.history'
2196 csa_bank=prefix(:lenpre)//'.CSA.bank'
2197 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2198 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2199 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2200 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2201 csa_int=prefix(:lenpre)//'.int'
2202 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2203 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2204 csa_in=prefix(:lenpre)//'.CSA.in'
2205 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2208 write (iout,'(80(1h-))')
2209 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2210 write (iout,'(80(1h-))')
2211 write (iout,*) "Input file : ",
2212 & pref_orig(:ilen(pref_orig))//'.inp'
2213 write (iout,*) "Output file : ",
2214 & outname(:ilen(outname))
2216 write (iout,*) "Sidechain potential file : ",
2217 & sidename(:ilen(sidename))
2219 write (iout,*) "SCp potential file : ",
2220 & scpname(:ilen(scpname))
2222 write (iout,*) "Electrostatic potential file : ",
2223 & elename(:ilen(elename))
2224 write (iout,*) "Cumulant coefficient file : ",
2225 & fouriername(:ilen(fouriername))
2226 write (iout,*) "Torsional parameter file : ",
2227 & torname(:ilen(torname))
2228 write (iout,*) "Double torsional parameter file : ",
2229 & tordname(:ilen(tordname))
2230 write (iout,*) "SCCOR parameter file : ",
2231 & sccorname(:ilen(sccorname))
2232 write (iout,*) "Bond & inertia constant file : ",
2233 & bondname(:ilen(bondname))
2234 write (iout,*) "Bending parameter file : ",
2235 & thetname(:ilen(thetname))
2236 write (iout,*) "Rotamer parameter file : ",
2237 & rotname(:ilen(rotname))
2238 write (iout,*) "Thetpdb parameter file : ",
2239 & thetname_pdb(:ilen(thetname_pdb))
2240 write (iout,*) "Threading database : ",
2241 & patname(:ilen(patname))
2243 &write (iout,*)" DIRTMP : ",
2245 write (iout,'(80(1h-))')
2249 c----------------------------------------------------------------------------
2250 subroutine card_concat(card)
2251 implicit real*8 (a-h,o-z)
2252 include 'DIMENSIONS'
2253 include 'COMMON.IOUNITS'
2255 character*80 karta,ucase
2257 read (inp,'(a)') karta
2260 do while (karta(80:80).eq.'&')
2261 card=card(:ilen(card)+1)//karta(:79)
2262 read (inp,'(a)') karta
2265 card=card(:ilen(card)+1)//karta
2268 c----------------------------------------------------------------------------------
2270 implicit real*8 (a-h,o-z)
2271 include 'DIMENSIONS'
2272 include 'COMMON.CHAIN'
2273 include 'COMMON.IOUNITS'
2275 open(irest2,file=rest2name,status='unknown')
2276 read(irest2,*) totT,EK,potE,totE,t_bath
2279 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2282 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2285 read (irest2,*) iset
2290 c---------------------------------------------------------------------------------
2291 subroutine read_fragments
2292 implicit real*8 (a-h,o-z)
2293 include 'DIMENSIONS'
2297 include 'COMMON.SETUP'
2298 include 'COMMON.CHAIN'
2299 include 'COMMON.IOUNITS'
2301 include 'COMMON.CONTROL'
2302 read(inp,*) nset,nfrag,npair,nfrag_back
2303 if(me.eq.king.or..not.out1file)
2304 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2305 & " nfrag_back",nfrag_back
2307 read(inp,*) mset(iset)
2309 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2311 if(me.eq.king.or..not.out1file)
2312 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2313 & ifrag(2,i,iset), qinfrag(i,iset)
2316 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2318 if(me.eq.king.or..not.out1file)
2319 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2320 & ipair(2,i,iset), qinpair(i,iset)
2323 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2324 & wfrag_back(3,i,iset),
2325 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2326 if(me.eq.king.or..not.out1file)
2327 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2328 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2333 C---------------------------------------------------------------------------
2334 subroutine read_afminp
2335 implicit real*8 (a-h,o-z)
2336 include 'DIMENSIONS'
2340 include 'COMMON.SETUP'
2341 include 'COMMON.CONTROL'
2342 include 'COMMON.CHAIN'
2343 include 'COMMON.IOUNITS'
2344 include 'COMMON.SBRIDGE'
2345 character*320 afmcard
2347 call card_concat(afmcard)
2348 call readi(afmcard,"BEG",afmbeg,0)
2349 call readi(afmcard,"END",afmend,0)
2350 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2351 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2352 print *,'FORCE=' ,forceAFMconst
2353 CCCC NOW PROPERTIES FOR AFM
2356 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2358 distafminit=dsqrt(distafminit)
2359 print *,'initdist',distafminit
2363 c-------------------------------------------------------------------------------
2364 subroutine read_dist_constr
2365 implicit real*8 (a-h,o-z)
2366 include 'DIMENSIONS'
2370 include 'COMMON.SETUP'
2371 include 'COMMON.CONTROL'
2372 include 'COMMON.CHAIN'
2373 include 'COMMON.IOUNITS'
2374 include 'COMMON.SBRIDGE'
2375 integer ifrag_(2,100),ipair_(2,100)
2376 double precision wfrag_(100),wpair_(100)
2377 character*500 controlcard
2379 write (iout,*) "Calling read_dist_constr"
2380 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2382 call card_concat(controlcard)
2383 call readi(controlcard,"NFRAG",nfrag_,0)
2384 call readi(controlcard,"NPAIR",npair_,0)
2385 call readi(controlcard,"NDIST",ndist_,0)
2386 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2387 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2388 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2389 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2390 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2391 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2392 c write (iout,*) "IFRAG"
2394 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2396 c write (iout,*) "IPAIR"
2398 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2402 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2403 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2404 & ifrag_(2,i)=nstart_sup+nsup-1
2405 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2407 if (wfrag_(i).gt.0.0d0) then
2408 do j=ifrag_(1,i),ifrag_(2,i)-1
2409 do k=j+1,ifrag_(2,i)
2410 c write (iout,*) "j",j," k",k
2412 if (constr_dist.eq.1) then
2417 forcon(nhpb)=wfrag_(i)
2418 else if (constr_dist.eq.2) then
2419 if (ddjk.le.dist_cut) then
2424 forcon(nhpb)=wfrag_(i)
2431 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2434 if (.not.out1file .or. me.eq.king)
2435 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2436 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2438 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2439 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2446 if (wpair_(i).gt.0.0d0) then
2454 do j=ifrag_(1,ii),ifrag_(2,ii)
2455 do k=ifrag_(1,jj),ifrag_(2,jj)
2459 forcon(nhpb)=wpair_(i)
2460 dhpb(nhpb)=dist(j,k)
2462 if (.not.out1file .or. me.eq.king)
2463 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2464 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2466 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2467 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2475 if (constr_dist.eq.11) then
2476 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2477 & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
2478 fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2481 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2482 & ibecarb(i),forcon(nhpb+1)
2484 if (forcon(nhpb+1).gt.0.0d0) then
2486 if (ibecarb(i).gt.0) then
2487 ihpb(i)=ihpb(i)+nres
2488 jhpb(i)=jhpb(i)+nres
2490 if (dhpb(nhpb).eq.0.0d0)
2491 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2493 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2494 C if (forcon(nhpb+1).gt.0.0d0) then
2496 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2498 if (.not.out1file .or. me.eq.king)
2499 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2500 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2502 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2503 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2510 c-------------------------------------------------------------------------------
2512 subroutine flush(iu)
2517 subroutine flush(iu)
2522 c------------------------------------------------------------------------------
2523 subroutine copy_to_tmp(source)
2524 include "DIMENSIONS"
2525 include "COMMON.IOUNITS"
2526 character*(*) source
2527 character* 256 tmpfile
2531 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2532 inquire(file=tmpfile,exist=ex)
2534 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2535 & " to temporary directory..."
2536 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2537 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2541 c------------------------------------------------------------------------------
2542 subroutine move_from_tmp(source)
2543 include "DIMENSIONS"
2544 include "COMMON.IOUNITS"
2545 character*(*) source
2548 write (*,*) "Moving ",source(:ilen(source)),
2549 & " from temporary directory to working directory"
2550 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2551 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2554 c------------------------------------------------------------------------------
2555 subroutine random_init(seed)
2557 C Initialize random number generator
2559 implicit real*8 (a-h,o-z)
2560 include 'DIMENSIONS'
2563 logical OKRandom, prng_restart
2565 integer iseed_array(4)
2567 include 'COMMON.IOUNITS'
2568 include 'COMMON.TIME1'
2569 include 'COMMON.THREAD'
2570 include 'COMMON.SBRIDGE'
2571 include 'COMMON.CONTROL'
2572 include 'COMMON.MCM'
2573 include 'COMMON.MAP'
2574 include 'COMMON.HEADER'
2575 include 'COMMON.CSA'
2576 include 'COMMON.CHAIN'
2577 include 'COMMON.MUCA'
2579 include 'COMMON.FFIELD'
2580 include 'COMMON.SETUP'
2581 iseed=-dint(dabs(seed))
2582 if (iseed.eq.0) then
2583 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2584 & 'Random seed undefined. The program will stop.'
2585 write (*,'(/80(1h*)/20x,a/80(1h*))')
2586 & 'Random seed undefined. The program will stop.'
2588 call mpi_finalize(mpi_comm_world,ierr)
2590 stop 'Bad random seed.'
2593 if (fg_rank.eq.0) then
2597 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2598 OKRandom = prng_restart(me,iseed)
2601 tmp=65536.0d0**(4-i)
2602 iseed_array(i) = dint(seed/tmp)
2603 seed=seed-iseed_array(i)*tmp
2606 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2607 & (iseed_array(i),i=1,4)
2608 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2609 & (iseed_array(i),i=1,4)
2610 OKRandom = prng_restart(me,iseed_array)
2613 c r1 = prng_next(me)
2614 r1=ran_number(0.0D0,1.0D0)
2616 & write (iout,*) 'ran_num',r1
2617 if (r1.lt.0.0d0) OKRandom=.false.
2619 if (.not.OKRandom) then
2620 write (iout,*) 'PRNG IS NOT WORKING!!!'
2621 print *,'PRNG IS NOT WORKING!!!'
2624 call mpi_abort(mpi_comm_world,error_msg,ierr)
2627 write (iout,*) 'too many processors for parallel prng'
2628 write (*,*) 'too many processors for parallel prng'
2636 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)