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 include 'COMMON.SHIELD'
85 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
86 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
88 character*320 controlcard
93 read (INP,'(a)') titel
94 call card_concat(controlcard)
95 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
96 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
97 call reada(controlcard,'SEED',seed,0.0D0)
98 call random_init(seed)
99 C Set up the time limit (caution! The time must be input in minutes!)
100 read_cart=index(controlcard,'READ_CART').gt.0
101 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
102 C this variable with_theta_constr is the variable which allow to read and execute the
103 C constrains on theta angles WITH_THETA_CONSTR is the keyword
104 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
105 write (iout,*) "with_theta_constr ",with_theta_constr
106 call readi(controlcard,'SYM',symetr,1)
107 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
108 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
109 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
110 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
111 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
112 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
113 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
114 call reada(controlcard,'DRMS',drms,0.1D0)
115 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
116 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
117 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
118 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
119 write (iout,'(a,f10.1)')'DRMS = ',drms
120 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
121 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
123 call readi(controlcard,'NZ_START',nz_start,0)
124 call readi(controlcard,'NZ_END',nz_end,0)
125 call readi(controlcard,'IZ_SC',iz_sc,0)
127 safety = 60.0d0*safety
130 call reada(controlcard,"T_BATH",t_bath,300.0d0)
131 minim=(index(controlcard,'MINIMIZE').gt.0)
132 dccart=(index(controlcard,'CART').gt.0)
133 overlapsc=(index(controlcard,'OVERLAP').gt.0)
134 overlapsc=.not.overlapsc
135 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
136 searchsc=.not.searchsc
137 sideadd=(index(controlcard,'SIDEADD').gt.0)
138 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
139 outpdb=(index(controlcard,'PDBOUT').gt.0)
140 outmol2=(index(controlcard,'MOL2OUT').gt.0)
141 pdbref=(index(controlcard,'PDBREF').gt.0)
142 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
143 indpdb=index(controlcard,'PDBSTART')
144 extconf=(index(controlcard,'EXTCONF').gt.0)
145 AFMlog=(index(controlcard,'AFM'))
146 selfguide=(index(controlcard,'SELFGUIDE'))
147 print *,'AFMlog',AFMlog,selfguide,"KUPA"
148 call readi(controlcard,'IPRINT',iprint,0)
149 C SHIELD keyword sets if the shielding effect of side-chains is used
150 C 0 denots no shielding is used all peptide are equally despite the
151 C solvent accesible area
152 C 1 the newly introduced function
153 C 2 reseved for further possible developement
154 call readi(controlcard,'SHIELD',shield_mode,0)
155 C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
156 write(iout,*) "shield_mode",shield_mode
158 call readi(controlcard,'MAXGEN',maxgen,10000)
159 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
160 call readi(controlcard,"KDIAG",kdiag,0)
161 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
162 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
163 & write (iout,*) "RESCALE_MODE",rescale_mode
164 split_ene=index(controlcard,'SPLIT_ENE').gt.0
165 if (index(controlcard,'REGULAR').gt.0.0D0) then
166 call reada(controlcard,'WEIDIS',weidis,0.1D0)
170 if (index(controlcard,'CHECKGRAD').gt.0) then
172 if (index(controlcard,'CART').gt.0) then
174 elseif (index(controlcard,'CARINT').gt.0) then
179 elseif (index(controlcard,'THREAD').gt.0) then
181 call readi(controlcard,'THREAD',nthread,0)
182 if (nthread.gt.0) then
183 call reada(controlcard,'WEIDIS',weidis,0.1D0)
186 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
187 stop 'Error termination in Read_Control.'
189 else if (index(controlcard,'MCMA').gt.0) then
191 else if (index(controlcard,'MCEE').gt.0) then
193 else if (index(controlcard,'MULTCONF').gt.0) then
195 else if (index(controlcard,'MAP').gt.0) then
197 call readi(controlcard,'MAP',nmap,0)
198 else if (index(controlcard,'CSA').gt.0) then
200 crc else if (index(controlcard,'ZSCORE').gt.0) then
202 crc ZSCORE is rm from UNRES, modecalc=9 is available
205 cfcm else if (index(controlcard,'MCMF').gt.0) then
207 else if (index(controlcard,'SOFTREG').gt.0) then
209 else if (index(controlcard,'CHECK_BOND').gt.0) then
211 else if (index(controlcard,'TEST').gt.0) then
213 else if (index(controlcard,'MD').gt.0) then
215 else if (index(controlcard,'RE ').gt.0) then
219 lmuca=index(controlcard,'MUCA').gt.0
220 call readi(controlcard,'MUCADYN',mucadyn,0)
221 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
222 if (lmuca .and. (me.eq.king .or. .not.out1file ))
224 write (iout,*) 'MUCADYN=',mucadyn
225 write (iout,*) 'MUCASMOOTH=',muca_smooth
228 iscode=index(controlcard,'ONE_LETTER')
229 indphi=index(controlcard,'PHI')
230 indback=index(controlcard,'BACK')
231 iranconf=index(controlcard,'RAND_CONF')
232 i2ndstr=index(controlcard,'USE_SEC_PRED')
233 gradout=index(controlcard,'GRADOUT').gt.0
234 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
235 C DISTCHAINMAX become obsolete for periodic boundry condition
236 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
237 C Reading the dimensions of box in x,y,z coordinates
238 call reada(controlcard,'BOXX',boxxsize,100.0d0)
239 call reada(controlcard,'BOXY',boxysize,100.0d0)
240 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
241 c Cutoff range for interactions
242 call reada(controlcard,"R_CUT",r_cut,15.0d0)
243 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
244 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
245 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
246 if (lipthick.gt.0.0d0) then
247 bordliptop=(boxzsize+lipthick)/2.0
248 bordlipbot=bordliptop-lipthick
250 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
251 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
252 buflipbot=bordlipbot+lipbufthick
253 bufliptop=bordliptop-lipbufthick
254 if ((lipbufthick*2.0d0).gt.lipthick)
255 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
257 write(iout,*) "bordliptop=",bordliptop
258 write(iout,*) "bordlipbot=",bordlipbot
259 write(iout,*) "bufliptop=",bufliptop
260 write(iout,*) "buflipbot=",buflipbot
261 write (iout,*) "SHIELD MODE",shield_mode
262 if (shield_mode.gt.0) then
264 C VSolvSphere the volume of solving sphere
266 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
267 C there will be no distinction between proline peptide group and normal peptide
268 C group in case of shielding parameters
269 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
270 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
271 write (iout,*) VSolvSphere,VSolvSphere_div
272 C long axis of side chain
274 long_r_sidechain(i)=vbldsc0(1,i)
275 short_r_sidechain(i)=sigma0(i)
279 if (me.eq.king .or. .not.out1file )
280 & write (iout,*) "DISTCHAINMAX",distchainmax
282 if(me.eq.king.or..not.out1file)
283 & write (iout,'(2a)') diagmeth(kdiag),
284 & ' routine used to diagonalize matrices.'
287 c--------------------------------------------------------------------------
288 subroutine read_REMDpar
292 implicit real*8 (a-h,o-z)
294 include 'COMMON.IOUNITS'
295 include 'COMMON.TIME1'
298 include 'COMMON.LANGEVIN'
300 include 'COMMON.LANGEVIN.lang0'
302 include 'COMMON.INTERACT'
303 include 'COMMON.NAMES'
305 include 'COMMON.REMD'
306 include 'COMMON.CONTROL'
307 include 'COMMON.SETUP'
309 character*320 controlcard
310 character*3200 controlcard1
311 integer iremd_m_total
313 if(me.eq.king.or..not.out1file)
314 & write (iout,*) "REMD setup"
316 call card_concat(controlcard)
317 call readi(controlcard,"NREP",nrep,3)
318 call readi(controlcard,"NSTEX",nstex,1000)
319 call reada(controlcard,"RETMIN",retmin,10.0d0)
320 call reada(controlcard,"RETMAX",retmax,1000.0d0)
321 mremdsync=(index(controlcard,'SYNC').gt.0)
322 call readi(controlcard,"NSYN",i_sync_step,100)
323 restart1file=(index(controlcard,'REST1FILE').gt.0)
324 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
325 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
326 if(max_cache_traj_use.gt.max_cache_traj)
327 & max_cache_traj_use=max_cache_traj
328 if(me.eq.king.or..not.out1file) then
329 cd if (traj1file) then
330 crc caching is in testing - NTWX is not ignored
331 cd write (iout,*) "NTWX value is ignored"
332 cd write (iout,*) " trajectory is stored to one file by master"
333 cd write (iout,*) " before exchange at NSTEX intervals"
335 write (iout,*) "NREP= ",nrep
336 write (iout,*) "NSTEX= ",nstex
337 write (iout,*) "SYNC= ",mremdsync
338 write (iout,*) "NSYN= ",i_sync_step
339 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
342 if (index(controlcard,'TLIST').gt.0) then
344 call card_concat(controlcard1)
345 read(controlcard1,*) (remd_t(i),i=1,nrep)
346 if(me.eq.king.or..not.out1file)
347 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
350 if (index(controlcard,'MLIST').gt.0) then
352 call card_concat(controlcard1)
353 read(controlcard1,*) (remd_m(i),i=1,nrep)
354 if(me.eq.king.or..not.out1file) then
355 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
358 iremd_m_total=iremd_m_total+remd_m(i)
360 write (iout,*) 'Total number of replicas ',iremd_m_total
363 if(me.eq.king.or..not.out1file)
364 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
367 c--------------------------------------------------------------------------
368 subroutine read_MDpar
372 implicit real*8 (a-h,o-z)
374 include 'COMMON.IOUNITS'
375 include 'COMMON.TIME1'
378 include 'COMMON.LANGEVIN'
380 include 'COMMON.LANGEVIN.lang0'
382 include 'COMMON.INTERACT'
383 include 'COMMON.NAMES'
385 include 'COMMON.SETUP'
386 include 'COMMON.CONTROL'
387 include 'COMMON.SPLITELE'
389 character*320 controlcard
391 call card_concat(controlcard)
392 call readi(controlcard,"NSTEP",n_timestep,1000000)
393 call readi(controlcard,"NTWE",ntwe,100)
394 call readi(controlcard,"NTWX",ntwx,1000)
395 call reada(controlcard,"DT",d_time,1.0d-1)
396 call reada(controlcard,"DVMAX",dvmax,2.0d1)
397 call reada(controlcard,"DAMAX",damax,1.0d1)
398 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
399 call readi(controlcard,"LANG",lang,0)
400 RESPA = index(controlcard,"RESPA") .gt. 0
401 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
402 ntime_split0=ntime_split
403 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
404 ntime_split0=ntime_split
405 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
406 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
407 rest = index(controlcard,"REST").gt.0
408 tbf = index(controlcard,"TBF").gt.0
409 usampl = index(controlcard,"USAMPL").gt.0
411 mdpdb = index(controlcard,"MDPDB").gt.0
412 call reada(controlcard,"T_BATH",t_bath,300.0d0)
413 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
414 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
415 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
416 if (count_reset_moment.eq.0) count_reset_moment=1000000000
417 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
418 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
419 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
420 if (count_reset_vel.eq.0) count_reset_vel=1000000000
421 large = index(controlcard,"LARGE").gt.0
422 print_compon = index(controlcard,"PRINT_COMPON").gt.0
423 rattle = index(controlcard,"RATTLE").gt.0
424 c if performing umbrella sampling, fragments constrained are read from the fragment file
430 if(me.eq.king.or..not.out1file) then
432 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
434 write (iout,'(a)') "The units are:"
435 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
436 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
437 & " acceleration: angstrom/(48.9 fs)**2"
438 write (iout,'(a)') "energy: kcal/mol, temperature: K"
440 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
441 write (iout,'(a60,f10.5,a)')
442 & "Initial time step of numerical integration:",d_time,
444 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
446 write (iout,'(2a,i4,a)')
447 & "A-MTS algorithm used; initial time step for fast-varying",
448 & " short-range forces split into",ntime_split," steps."
449 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
450 & r_cut," lambda",rlamb
452 write (iout,'(2a,f10.5)')
453 & "Maximum acceleration threshold to reduce the time step",
454 & "/increase split number:",damax
455 write (iout,'(2a,f10.5)')
456 & "Maximum predicted energy drift to reduce the timestep",
457 & "/increase split number:",edriftmax
458 write (iout,'(a60,f10.5)')
459 & "Maximum velocity threshold to reduce velocities:",dvmax
460 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
461 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
462 if (rattle) write (iout,'(a60)')
463 & "Rattle algorithm used to constrain the virtual bonds"
467 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
468 call reada(controlcard,"RWAT",rwat,1.4d0)
469 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
470 surfarea=index(controlcard,"SURFAREA").gt.0
471 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
472 if(me.eq.king.or..not.out1file)then
473 write (iout,'(/a,$)') "Langevin dynamics calculation"
476 & " with direct integration of Langevin equations"
477 else if (lang.eq.2) then
478 write (iout,'(a/)') " with TINKER stochasic MD integrator"
479 else if (lang.eq.3) then
480 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
481 else if (lang.eq.4) then
482 write (iout,'(a/)') " in overdamped mode"
484 write (iout,'(//a,i5)')
485 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
488 write (iout,'(a60,f10.5)') "Temperature:",t_bath
489 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
490 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
491 write (iout,'(a60,f10.5)')
492 & "Scaling factor of the friction forces:",scal_fric
493 if (surfarea) write (iout,'(2a,i10,a)')
494 & "Friction coefficients will be scaled by solvent-accessible",
495 & " surface area every",reset_fricmat," steps."
497 c Calculate friction coefficients and bounds of stochastic forces
498 eta=6*pi*cPoise*etawat
499 if(me.eq.king.or..not.out1file)
500 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
502 gamp=scal_fric*(pstok+rwat)*eta
503 stdfp=dsqrt(2*Rb*t_bath/d_time)
505 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
506 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
508 if(me.eq.king.or..not.out1file)then
509 write (iout,'(/2a/)')
510 & "Radii of site types and friction coefficients and std's of",
511 & " stochastic forces of fully exposed sites"
512 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
514 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
515 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
519 if(me.eq.king.or..not.out1file)then
520 write (iout,'(a)') "Berendsen bath calculation"
521 write (iout,'(a60,f10.5)') "Temperature:",t_bath
522 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
524 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
525 & count_reset_moment," steps"
527 & write (iout,'(a,i10,a)')
528 & "Velocities will be reset at random every",count_reset_vel,
532 if(me.eq.king.or..not.out1file)
533 & write (iout,'(a31)') "Microcanonical mode calculation"
535 if(me.eq.king.or..not.out1file)then
536 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
538 write(iout,*) "MD running with constraints."
539 write(iout,*) "Equilibration time ", eq_time, " mtus."
540 write(iout,*) "Constraining ", nfrag," fragments."
541 write(iout,*) "Length of each fragment, weight and q0:"
543 write (iout,*) "Set of restraints #",iset
545 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
546 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
548 write(iout,*) "constraints between ", npair, "fragments."
549 write(iout,*) "constraint pairs, weights and q0:"
551 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
552 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
554 write(iout,*) "angle constraints within ", nfrag_back,
555 & "backbone fragments."
556 write(iout,*) "fragment, weights:"
558 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
559 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
560 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
563 iset=mod(kolor,nset)+1
566 if(me.eq.king.or..not.out1file)
567 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
570 c------------------------------------------------------------------------------
573 C Read molecular data.
575 implicit real*8 (a-h,o-z)
581 include 'COMMON.IOUNITS'
584 include 'COMMON.INTERACT'
585 include 'COMMON.LOCAL'
586 include 'COMMON.NAMES'
587 include 'COMMON.CHAIN'
588 include 'COMMON.FFIELD'
589 include 'COMMON.SBRIDGE'
590 include 'COMMON.HEADER'
591 include 'COMMON.CONTROL'
592 include 'COMMON.DBASE'
593 include 'COMMON.THREAD'
594 include 'COMMON.CONTACTS'
595 include 'COMMON.TORCNSTR'
596 include 'COMMON.TIME1'
597 include 'COMMON.BOUNDS'
599 include 'COMMON.SETUP'
600 character*4 sequence(maxres)
602 double precision x(maxvar)
603 character*256 pdbfile
604 character*400 weightcard
605 character*80 weightcard_t,ucase
606 dimension itype_pdb(maxres)
607 common /pizda/ itype_pdb
608 logical seq_comp,fail
609 double precision energia(0:n_ene)
615 C Read weights of the subsequent energy terms.
616 call card_concat(weightcard)
617 call reada(weightcard,'WLONG',wlong,1.0D0)
618 call reada(weightcard,'WSC',wsc,wlong)
619 call reada(weightcard,'WSCP',wscp,wlong)
620 call reada(weightcard,'WELEC',welec,1.0D0)
621 call reada(weightcard,'WVDWPP',wvdwpp,welec)
622 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
623 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
624 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
625 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
626 call reada(weightcard,'WTURN3',wturn3,1.0D0)
627 call reada(weightcard,'WTURN4',wturn4,1.0D0)
628 call reada(weightcard,'WTURN6',wturn6,1.0D0)
629 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
630 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
631 call reada(weightcard,'WBOND',wbond,1.0D0)
632 call reada(weightcard,'WTOR',wtor,1.0D0)
633 call reada(weightcard,'WTORD',wtor_d,1.0D0)
634 call reada(weightcard,'WANG',wang,1.0D0)
635 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
636 call reada(weightcard,'SCAL14',scal14,0.4D0)
637 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
638 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
639 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
640 call reada(weightcard,'TEMP0',temp0,300.0d0)
641 call reada(weightcard,'WLT',wliptran,0.0D0)
642 if (index(weightcard,'SOFT').gt.0) ipot=6
643 C 12/1/95 Added weight for the multi-body term WCORR
644 call reada(weightcard,'WCORRH',wcorr,1.0D0)
645 if (wcorr4.gt.0.0d0) wcorr=wcorr4
665 if(me.eq.king.or..not.out1file)
666 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
667 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
669 10 format (/'Energy-term weights (unscaled):'//
670 & 'WSCC= ',f10.6,' (SC-SC)'/
671 & 'WSCP= ',f10.6,' (SC-p)'/
672 & 'WELEC= ',f10.6,' (p-p electr)'/
673 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
674 & 'WBOND= ',f10.6,' (stretching)'/
675 & 'WANG= ',f10.6,' (bending)'/
676 & 'WSCLOC= ',f10.6,' (SC local)'/
677 & 'WTOR= ',f10.6,' (torsional)'/
678 & 'WTORD= ',f10.6,' (double torsional)'/
679 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
680 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
681 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
682 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
683 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
684 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
685 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
686 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
687 & 'WTURN6= ',f10.6,' (turns, 6th order)')
688 if(me.eq.king.or..not.out1file)then
689 if (wcorr4.gt.0.0d0) then
690 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
691 & 'between contact pairs of peptide groups'
692 write (iout,'(2(a,f5.3/))')
693 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
694 & 'Range of quenching the correlation terms:',2*delt_corr
695 else if (wcorr.gt.0.0d0) then
696 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
697 & 'between contact pairs of peptide groups'
699 write (iout,'(a,f8.3)')
700 & 'Scaling factor of 1,4 SC-p interactions:',scal14
701 write (iout,'(a,f8.3)')
702 & 'General scaling factor of SC-p interactions:',scalscp
704 r0_corr=cutoff_corr-delt_corr
706 aad(i,1)=scalscp*aad(i,1)
707 aad(i,2)=scalscp*aad(i,2)
708 bad(i,1)=scalscp*bad(i,1)
709 bad(i,2)=scalscp*bad(i,2)
711 call rescale_weights(t_bath)
712 if(me.eq.king.or..not.out1file)
713 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
714 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
716 22 format (/'Energy-term weights (scaled):'//
717 & 'WSCC= ',f10.6,' (SC-SC)'/
718 & 'WSCP= ',f10.6,' (SC-p)'/
719 & 'WELEC= ',f10.6,' (p-p electr)'/
720 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
721 & 'WBOND= ',f10.6,' (stretching)'/
722 & 'WANG= ',f10.6,' (bending)'/
723 & 'WSCLOC= ',f10.6,' (SC local)'/
724 & 'WTOR= ',f10.6,' (torsional)'/
725 & 'WTORD= ',f10.6,' (double torsional)'/
726 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
727 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
728 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
729 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
730 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
731 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
732 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
733 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
734 & 'WTURN6= ',f10.6,' (turns, 6th order)')
735 if(me.eq.king.or..not.out1file)
736 & write (iout,*) "Reference temperature for weights calculation:",
738 call reada(weightcard,"D0CM",d0cm,3.78d0)
739 call reada(weightcard,"AKCM",akcm,15.1d0)
740 call reada(weightcard,"AKTH",akth,11.0d0)
741 call reada(weightcard,"AKCT",akct,12.0d0)
742 call reada(weightcard,"V1SS",v1ss,-1.08d0)
743 call reada(weightcard,"V2SS",v2ss,7.61d0)
744 call reada(weightcard,"V3SS",v3ss,13.7d0)
745 call reada(weightcard,"EBR",ebr,-5.50D0)
746 call reada(weightcard,"ATRISS",atriss,0.301D0)
747 call reada(weightcard,"BTRISS",btriss,0.021D0)
748 call reada(weightcard,"CTRISS",ctriss,1.001D0)
749 call reada(weightcard,"DTRISS",dtriss,1.001D0)
750 write (iout,*) "ATRISS=", atriss
751 write (iout,*) "BTRISS=", btriss
752 write (iout,*) "CTRISS=", ctriss
753 write (iout,*) "DTRISS=", dtriss
754 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
756 dyn_ss_mask(i)=.false.
760 dyn_ssbond_ij(i,j)=1.0d300
763 call reada(weightcard,"HT",Ht,0.0D0)
765 ss_depth=ebr/wsc-0.25*eps(1,1)
766 Ht=Ht/wsc-0.25*eps(1,1)
767 akcm=akcm*wstrain/wsc
768 akth=akth*wstrain/wsc
769 akct=akct*wstrain/wsc
770 v1ss=v1ss*wstrain/wsc
771 v2ss=v2ss*wstrain/wsc
772 v3ss=v3ss*wstrain/wsc
774 if (wstrain.ne.0.0) then
775 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
781 if(me.eq.king.or..not.out1file) then
782 write (iout,*) "Parameters of the SS-bond potential:"
783 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
785 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
786 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
787 write (iout,*)" HT",Ht
788 print *,'indpdb=',indpdb,' pdbref=',pdbref
790 if (indpdb.gt.0 .or. pdbref) then
791 read(inp,'(a)') pdbfile
792 if(me.eq.king.or..not.out1file)
793 & write (iout,'(2a)') 'PDB data will be read from file ',
794 & pdbfile(:ilen(pdbfile))
795 open(ipdbin,file=pdbfile,status='old',err=33)
797 33 write (iout,'(a)') 'Error opening PDB file.'
800 c write (iout,*) 'Begin reading pdb data'
803 c write (iout,*) 'Finished reading pdb data'
805 if(me.eq.king.or..not.out1file)
806 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
807 & ' nstart_sup=',nstart_sup
809 itype_pdb(i)=itype(i)
813 nct=nstart_sup+nsup-1
814 call contact(.false.,ncont_ref,icont_ref,co)
817 if(me.eq.king.or..not.out1file)
818 & write(iout,*)'Adding sidechains'
822 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
825 do while (fail.and.nsi.le.maxsi)
826 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
829 if(fail) write(iout,*)'Adding sidechain failed for res ',
830 & i,' after ',nsi,' trials'
835 if (indpdb.eq.0) then
836 C Read sequence if not taken from the pdb file.
838 c print *,'nres=',nres
839 if (iscode.gt.0) then
840 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
842 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
844 C Convert sequence to numeric code
846 itype(i)=rescode(i,sequence(i),iscode)
848 C Assign initial virtual bond lengths
854 vbld(i+nres)=dsc(iabs(itype(i)))
855 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
856 c write (iout,*) "i",i," itype",itype(i),
857 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
861 c print '(20i4)',(itype(i),i=1,nres)
864 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
866 if (itype(i).eq.ntyp1) then
870 else if (iabs(itype(i+1)).ne.20) then
872 else if (iabs(itype(i)).ne.20) then
879 if(me.eq.king.or..not.out1file)then
880 write (iout,*) "ITEL"
882 write (iout,*) i,itype(i),itel(i)
884 print *,'Call Read_Bridge.'
887 C 8/13/98 Set limits to generating the dihedral angles
892 read (inp,*) ndih_constr
893 if (ndih_constr.gt.0) then
895 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
897 if(me.eq.king.or..not.out1file)then
899 & 'There are',ndih_constr,' constraints on phi angles.'
901 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
906 phi0(i)=deg2rad*phi0(i)
907 drange(i)=deg2rad*drange(i)
909 C if(me.eq.king.or..not.out1file)
910 C & write (iout,*) 'FTORS',ftors
913 phibound(1,ii) = phi0(i)-drange(i)
914 phibound(2,ii) = phi0(i)+drange(i)
917 C first setting the theta boundaries to 0 to pi
918 C this mean that there is no energy penalty for any angle occuring this can be applied
919 C for generate random conformation but is not implemented in this way
924 C begin reading theta constrains this is quartic constrains allowing to
925 C have smooth second derivative
926 if (with_theta_constr) then
927 C with_theta_constr is keyword allowing for occurance of theta constrains
928 read (inp,*) ntheta_constr
929 C ntheta_constr is the number of theta constrains
930 if (ntheta_constr.gt.0) then
932 read (inp,*) (itheta_constr(i),theta_constr0(i),
933 & theta_drange(i),for_thet_constr(i),
935 C the above code reads from 1 to ntheta_constr
936 C itheta_constr(i) residue i for which is theta_constr
937 C theta_constr0 the global minimum value
938 C theta_drange is range for which there is no energy penalty
939 C for_thet_constr is the force constant for quartic energy penalty
941 if(me.eq.king.or..not.out1file)then
943 & 'There are',ntheta_constr,' constraints on phi angles.'
945 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
951 theta_constr0(i)=deg2rad*theta_constr0(i)
952 theta_drange(i)=deg2rad*theta_drange(i)
954 C if(me.eq.king.or..not.out1file)
955 C & write (iout,*) 'FTORS',ftors
956 C do i=1,ntheta_constr
957 C ii = itheta_constr(i)
958 C thetabound(1,ii) = phi0(i)-drange(i)
959 C thetabound(2,ii) = phi0(i)+drange(i)
961 endif ! ntheta_constr.gt.0
962 endif! with_theta_constr
964 C with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
965 C write (iout,*) "with_dihed_constr ",with_dihed_constr
970 write (iout,'(a)') 'Boundaries in phi angle sampling:'
972 write (iout,'(a3,i5,2f10.1)')
973 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
979 cd print *,'NNT=',NNT,' NCT=',NCT
980 if (itype(1).eq.ntyp1) nnt=2
981 if (itype(nres).eq.ntyp1) nct=nct-1
983 if(me.eq.king.or..not.out1file)
984 & write (iout,'(a,i3)') 'nsup=',nsup
986 if (nsup.le.(nct-nnt+1)) then
987 do i=0,nct-nnt+1-nsup
988 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
994 & 'Error - sequences to be superposed do not match.'
997 do i=0,nsup-(nct-nnt+1)
998 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1000 nstart_sup=nstart_sup+i
1006 & 'Error - sequences to be superposed do not match.'
1009 if (nsup.eq.0) nsup=nct-nnt
1010 if (nstart_sup.eq.0) nstart_sup=nnt
1011 if (nstart_seq.eq.0) nstart_seq=nnt
1012 if(me.eq.king.or..not.out1file)
1013 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1014 & ' nstart_seq=',nstart_seq
1016 c--- Zscore rms -------
1017 if (nz_start.eq.0) nz_start=nnt
1018 if (nz_end.eq.0 .and. nsup.gt.0) then
1020 else if (nz_end.eq.0) then
1023 if(me.eq.king.or..not.out1file)then
1024 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1025 write (iout,*) 'IZ_SC=',iz_sc
1027 c----------------------
1030 if (.not.pdbref) then
1031 call read_angles(inp,*38)
1033 38 write (iout,'(a)') 'Error reading reference structure.'
1035 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1036 stop 'Error reading reference structure'
1040 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1047 cref(j,i,kkk)=c(j,i)
1050 call contact(.true.,ncont_ref,icont_ref,co)
1054 write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1056 if (constr_dist.gt.0) call read_dist_constr
1057 write (iout,*) "After read_dist_constr nhpb",nhpb
1058 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1060 if(me.eq.king.or..not.out1file)
1061 & write (iout,*) 'Contact order:',co
1063 if(me.eq.king.or..not.out1file)
1064 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1067 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1069 if(me.eq.king.or..not.out1file)
1070 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1071 & icont_ref(1,i),' ',
1072 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1076 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1077 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1078 & modecalc.ne.10) then
1079 C If input structure hasn't been supplied from the PDB file read or generate
1081 if (iranconf.eq.0 .and. .not. extconf) then
1082 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1083 & write (iout,'(a)') 'Initial geometry will be read in.'
1085 read(inp,'(8f10.5)',end=36,err=36)
1086 & ((c(l,k),l=1,3),k=1,nres),
1087 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1088 write (iout,*) "Exit READ_CART"
1089 write (iout,'(8f10.5)')
1090 & ((c(l,k),l=1,3),k=1,nres),
1091 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1092 call int_from_cart1(.true.)
1093 write (iout,*) "Finish INT_TO_CART"
1096 dc(j,i)=c(j,i+1)-c(j,i)
1097 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1101 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1103 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1104 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1110 call read_angles(inp,*36)
1113 36 write (iout,'(a)') 'Error reading angle file.'
1115 call mpi_finalize( MPI_COMM_WORLD,IERR )
1117 stop 'Error reading angle file.'
1119 else if (extconf) then
1120 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1121 & write (iout,'(a)') 'Extended chain initial geometry.'
1123 theta(i)=90d0*deg2rad
1126 phi(i)=180d0*deg2rad
1129 alph(i)=110d0*deg2rad
1132 omeg(i)=-120d0*deg2rad
1133 if (itype(i).le.0) omeg(i)=-omeg(i)
1135 call chainbuild_extconf
1137 if(me.eq.king.or..not.out1file)
1138 & write (iout,'(a)') 'Random-generated initial geometry.'
1142 if (me.eq.king .or. fg_rank.eq.0 .and. (
1143 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1147 call gen_rand_conf(itmp,*30)
1149 30 write (iout,*) 'Failed to generate random conformation',
1150 & ', itrial=',itrial
1151 write (*,*) 'Processor:',me,
1152 & ' Failed to generate random conformation',
1162 write (iout,'(a,i3,a)') 'Processor:',me,
1163 & ' error in generating random conformation.'
1164 write (*,'(a,i3,a)') 'Processor:',me,
1165 & ' error in generating random conformation.'
1168 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1173 & ' error in generating random conformation.'
1178 elseif (modecalc.eq.4) then
1179 read (inp,'(a)') intinname
1180 open (intin,file=intinname,status='old',err=333)
1181 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1182 & write (iout,'(a)') 'intinname',intinname
1183 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1185 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1187 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1189 stop 'Error opening angle file.'
1193 C Generate distance constraints, if the PDB structure is to be regularized.
1194 if (nthread.gt.0) then
1195 call read_threadbase
1198 if (me.eq.king .or. .not. out1file)
1200 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1201 write (iout,'(/a,i3,a)')
1202 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1203 write (iout,'(20i4)') (iss(i),i=1,ns)
1205 write(iout,*)"Running with dynamic disulfide-bond formation"
1207 write (iout,'(/a/)') 'Pre-formed links are:'
1213 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1214 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1220 if (ns.gt.0.and.dyn_ss) then
1224 forcon(i-nss)=forcon(i)
1231 dyn_ss_mask(iss(i))=.true.
1234 if (i2ndstr.gt.0) call secstrp2dihc
1235 c call geom_to_var(nvar,x)
1236 c call etotal(energia(0))
1237 c call enerprint(energia(0))
1238 c call briefout(0,etot)
1240 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1241 cd write (iout,'(a)') 'Variable list:'
1242 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1244 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1245 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1246 & 'Processor',myrank,': end reading molecular data.'
1251 c--------------------------------------------------------------------------
1252 logical function seq_comp(itypea,itypeb,length)
1254 integer length,itypea(length),itypeb(length)
1257 if (itypea(i).ne.itypeb(i)) then
1265 c-----------------------------------------------------------------------------
1266 subroutine read_bridge
1267 C Read information about disulfide bridges.
1268 implicit real*8 (a-h,o-z)
1269 include 'DIMENSIONS'
1273 include 'COMMON.IOUNITS'
1274 include 'COMMON.GEO'
1275 include 'COMMON.VAR'
1276 include 'COMMON.INTERACT'
1277 include 'COMMON.LOCAL'
1278 include 'COMMON.NAMES'
1279 include 'COMMON.CHAIN'
1280 include 'COMMON.FFIELD'
1281 include 'COMMON.SBRIDGE'
1282 include 'COMMON.HEADER'
1283 include 'COMMON.CONTROL'
1284 include 'COMMON.DBASE'
1285 include 'COMMON.THREAD'
1286 include 'COMMON.TIME1'
1287 include 'COMMON.SETUP'
1288 C Read bridging residues.
1289 read (inp,*) ns,(iss(i),i=1,ns)
1291 if(me.eq.king.or..not.out1file)
1292 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1293 C Check whether the specified bridging residues are cystines.
1295 if (itype(iss(i)).ne.1) then
1296 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1297 & 'Do you REALLY think that the residue ',
1298 & restyp(itype(iss(i))),i,
1299 & ' can form a disulfide bridge?!!!'
1300 write (*,'(2a,i3,a)')
1301 & 'Do you REALLY think that the residue ',
1302 & restyp(itype(iss(i))),i,
1303 & ' can form a disulfide bridge?!!!'
1305 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1310 C Read preformed bridges.
1312 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1314 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1317 C Check if the residues involved in bridges are in the specified list of
1318 C bridging residues.
1321 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1322 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1323 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1324 & ' contains residues present in other pairs.'
1325 write (*,'(a,i3,a)') 'Disulfide pair',i,
1326 & ' contains residues present in other pairs.'
1328 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1334 if (ihpb(i).eq.iss(j)) goto 10
1336 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1339 if (jhpb(i).eq.iss(j)) goto 20
1341 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1347 ihpb(i)=ihpb(i)+nres
1348 jhpb(i)=jhpb(i)+nres
1354 c----------------------------------------------------------------------------
1355 subroutine read_x(kanal,*)
1356 implicit real*8 (a-h,o-z)
1357 include 'DIMENSIONS'
1358 include 'COMMON.GEO'
1359 include 'COMMON.VAR'
1360 include 'COMMON.CHAIN'
1361 include 'COMMON.IOUNITS'
1362 include 'COMMON.CONTROL'
1363 include 'COMMON.LOCAL'
1364 include 'COMMON.INTERACT'
1365 c Read coordinates from input
1367 read(kanal,'(8f10.5)',end=10,err=10)
1368 & ((c(l,k),l=1,3),k=1,nres),
1369 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1372 c(j,2*nres)=c(j,nres)
1374 call int_from_cart1(.false.)
1377 dc(j,i)=c(j,i+1)-c(j,i)
1378 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1382 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1384 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1385 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1393 c----------------------------------------------------------------------------
1394 subroutine read_threadbase
1395 implicit real*8 (a-h,o-z)
1396 include 'DIMENSIONS'
1397 include 'COMMON.IOUNITS'
1398 include 'COMMON.GEO'
1399 include 'COMMON.VAR'
1400 include 'COMMON.INTERACT'
1401 include 'COMMON.LOCAL'
1402 include 'COMMON.NAMES'
1403 include 'COMMON.CHAIN'
1404 include 'COMMON.FFIELD'
1405 include 'COMMON.SBRIDGE'
1406 include 'COMMON.HEADER'
1407 include 'COMMON.CONTROL'
1408 include 'COMMON.DBASE'
1409 include 'COMMON.THREAD'
1410 include 'COMMON.TIME1'
1411 C Read pattern database for threading.
1412 read (icbase,*) nseq
1414 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1415 & nres_base(2,i),nres_base(3,i)
1416 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1418 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1419 c & nres_base(2,i),nres_base(3,i)
1420 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1424 if (weidis.eq.0.0D0) weidis=0.1D0
1433 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1434 write (iout,'(a,i5)') 'nexcl: ',nexcl
1435 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1438 c------------------------------------------------------------------------------
1439 subroutine setup_var
1440 implicit real*8 (a-h,o-z)
1441 include 'DIMENSIONS'
1442 include 'COMMON.IOUNITS'
1443 include 'COMMON.GEO'
1444 include 'COMMON.VAR'
1445 include 'COMMON.INTERACT'
1446 include 'COMMON.LOCAL'
1447 include 'COMMON.NAMES'
1448 include 'COMMON.CHAIN'
1449 include 'COMMON.FFIELD'
1450 include 'COMMON.SBRIDGE'
1451 include 'COMMON.HEADER'
1452 include 'COMMON.CONTROL'
1453 include 'COMMON.DBASE'
1454 include 'COMMON.THREAD'
1455 include 'COMMON.TIME1'
1456 C Set up variable list.
1462 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1464 ialph(i,1)=nvar+nside
1468 if (indphi.gt.0) then
1470 else if (indback.gt.0) then
1475 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1478 c----------------------------------------------------------------------------
1479 subroutine gen_dist_constr
1480 C Generate CA distance constraints.
1481 implicit real*8 (a-h,o-z)
1482 include 'DIMENSIONS'
1483 include 'COMMON.IOUNITS'
1484 include 'COMMON.GEO'
1485 include 'COMMON.VAR'
1486 include 'COMMON.INTERACT'
1487 include 'COMMON.LOCAL'
1488 include 'COMMON.NAMES'
1489 include 'COMMON.CHAIN'
1490 include 'COMMON.FFIELD'
1491 include 'COMMON.SBRIDGE'
1492 include 'COMMON.HEADER'
1493 include 'COMMON.CONTROL'
1494 include 'COMMON.DBASE'
1495 include 'COMMON.THREAD'
1496 include 'COMMON.TIME1'
1497 dimension itype_pdb(maxres)
1498 common /pizda/ itype_pdb
1500 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1501 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1502 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1504 do i=nstart_sup,nstart_sup+nsup-1
1505 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1506 cd & ' seq_pdb', restyp(itype_pdb(i))
1507 do j=i+2,nstart_sup+nsup-1
1509 ihpb(nhpb)=i+nstart_seq-nstart_sup
1510 jhpb(nhpb)=j+nstart_seq-nstart_sup
1512 dhpb(nhpb)=dist(i,j)
1515 cd write (iout,'(a)') 'Distance constraints:'
1520 cd if (ii.gt.nres) then
1525 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1526 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1527 cd & dhpb(i),forcon(i)
1531 c----------------------------------------------------------------------------
1533 implicit real*8 (a-h,o-z)
1534 include 'DIMENSIONS'
1535 include 'COMMON.MAP'
1536 include 'COMMON.IOUNITS'
1537 character*3 angid(4) /'THE','PHI','ALP','OME'/
1538 character*80 mapcard,ucase
1540 read (inp,'(a)') mapcard
1541 mapcard=ucase(mapcard)
1542 if (index(mapcard,'PHI').gt.0) then
1544 else if (index(mapcard,'THE').gt.0) then
1546 else if (index(mapcard,'ALP').gt.0) then
1548 else if (index(mapcard,'OME').gt.0) then
1551 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1552 stop 'Error - illegal variable spec in MAP card.'
1554 call readi (mapcard,'RES1',res1(imap),0)
1555 call readi (mapcard,'RES2',res2(imap),0)
1556 if (res1(imap).eq.0) then
1557 res1(imap)=res2(imap)
1558 else if (res2(imap).eq.0) then
1559 res2(imap)=res1(imap)
1561 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1563 & 'Error - illegal definition of variable group in MAP.'
1564 stop 'Error - illegal definition of variable group in MAP.'
1566 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1567 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1568 call readi(mapcard,'NSTEP',nstep(imap),0)
1569 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1571 & 'Illegal boundary and/or step size specification in MAP.'
1572 stop 'Illegal boundary and/or step size specification in MAP.'
1577 c----------------------------------------------------------------------------
1579 implicit real*8 (a-h,o-z)
1580 include 'DIMENSIONS'
1581 include 'COMMON.IOUNITS'
1582 include 'COMMON.GEO'
1583 include 'COMMON.CSA'
1584 include 'COMMON.BANK'
1585 include 'COMMON.CONTROL'
1587 character*620 mcmcard
1588 call card_concat(mcmcard)
1590 call readi(mcmcard,'NCONF',nconf,50)
1591 call readi(mcmcard,'NADD',nadd,0)
1592 call readi(mcmcard,'JSTART',jstart,1)
1593 call readi(mcmcard,'JEND',jend,1)
1594 call readi(mcmcard,'NSTMAX',nstmax,500000)
1595 call readi(mcmcard,'N0',n0,1)
1596 call readi(mcmcard,'N1',n1,6)
1597 call readi(mcmcard,'N2',n2,4)
1598 call readi(mcmcard,'N3',n3,0)
1599 call readi(mcmcard,'N4',n4,0)
1600 call readi(mcmcard,'N5',n5,0)
1601 call readi(mcmcard,'N6',n6,10)
1602 call readi(mcmcard,'N7',n7,0)
1603 call readi(mcmcard,'N8',n8,0)
1604 call readi(mcmcard,'N9',n9,0)
1605 call readi(mcmcard,'N14',n14,0)
1606 call readi(mcmcard,'N15',n15,0)
1607 call readi(mcmcard,'N16',n16,0)
1608 call readi(mcmcard,'N17',n17,0)
1609 call readi(mcmcard,'N18',n18,0)
1611 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1613 call readi(mcmcard,'NDIFF',ndiff,2)
1614 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1615 call readi(mcmcard,'IS1',is1,1)
1616 call readi(mcmcard,'IS2',is2,8)
1617 call readi(mcmcard,'NRAN0',nran0,4)
1618 call readi(mcmcard,'NRAN1',nran1,2)
1619 call readi(mcmcard,'IRR',irr,1)
1620 call readi(mcmcard,'NSEED',nseed,20)
1621 call readi(mcmcard,'NTOTAL',ntotal,10000)
1622 call reada(mcmcard,'CUT1',cut1,2.0d0)
1623 call reada(mcmcard,'CUT2',cut2,5.0d0)
1624 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1625 call readi(mcmcard,'ICMAX',icmax,3)
1626 call readi(mcmcard,'IRESTART',irestart,0)
1627 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1630 call reada(mcmcard,'DELE',dele,20.0d0)
1631 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1632 call readi(mcmcard,'IREF',iref,0)
1633 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1634 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1635 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1636 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1637 write (iout,*) "NCONF_IN",nconf_in
1640 c----------------------------------------------------------------------------
1641 cfmc subroutine mcmfread
1642 cfmc implicit real*8 (a-h,o-z)
1643 cfmc include 'DIMENSIONS'
1644 cfmc include 'COMMON.MCMF'
1645 cfmc include 'COMMON.IOUNITS'
1646 cfmc include 'COMMON.GEO'
1647 cfmc character*80 ucase
1648 cfmc character*620 mcmcard
1649 cfmc call card_concat(mcmcard)
1651 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1652 cfmc write(iout,*)'MAXRANT=',maxrant
1653 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1654 cfmc write(iout,*)'MAXFAM=',maxfam
1655 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1656 cfmc write(iout,*)'NNET1=',nnet1
1657 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1658 cfmc write(iout,*)'NNET2=',nnet2
1659 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1660 cfmc write(iout,*)'NNET3=',nnet3
1661 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1662 cfmc write(iout,*)'ILASTT=',ilastt
1663 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1664 cfmc write(iout,*)'MAXSTR=',maxstr
1665 cfmc maxstr_f=maxstr/maxfam
1666 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1667 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1668 cfmc write(iout,*)'NMCMF=',nmcmf
1669 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1670 cfmc write(iout,*)'IFOCUS=',ifocus
1671 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1672 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1673 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1674 cfmc write(iout,*)'INTPRT=',intprt
1675 cfmc call readi(mcmcard,'IPRT',iprt,100)
1676 cfmc write(iout,*)'IPRT=',iprt
1677 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1678 cfmc write(iout,*)'IMAXTR=',imaxtr
1679 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1680 cfmc write(iout,*)'MAXEVEN=',maxeven
1681 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1682 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1683 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1684 cfmc write(iout,*)'INIMIN=',inimin
1685 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1686 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1687 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1688 cfmc write(iout,*)'NTHREAD=',nthread
1689 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1690 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1691 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1692 cfmc write(iout,*)'MAXPERT=',maxpert
1693 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1694 cfmc write(iout,*)'IRMSD=',irmsd
1695 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1696 cfmc write(iout,*)'DENEMIN=',denemin
1697 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1698 cfmc write(iout,*)'RCUT1S=',rcut1s
1699 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1700 cfmc write(iout,*)'RCUT1E=',rcut1e
1701 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1702 cfmc write(iout,*)'RCUT2S=',rcut2s
1703 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1704 cfmc write(iout,*)'RCUT2E=',rcut2e
1705 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1706 cfmc write(iout,*)'DPERT1=',d_pert1
1707 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1708 cfmc write(iout,*)'DPERT1A=',d_pert1a
1709 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1710 cfmc write(iout,*)'DPERT2=',d_pert2
1711 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1712 cfmc write(iout,*)'DPERT2A=',d_pert2a
1713 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1714 cfmc write(iout,*)'DPERT2B=',d_pert2b
1715 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1716 cfmc write(iout,*)'DPERT2C=',d_pert2c
1717 cfmc d_pert1=deg2rad*d_pert1
1718 cfmc d_pert1a=deg2rad*d_pert1a
1719 cfmc d_pert2=deg2rad*d_pert2
1720 cfmc d_pert2a=deg2rad*d_pert2a
1721 cfmc d_pert2b=deg2rad*d_pert2b
1722 cfmc d_pert2c=deg2rad*d_pert2c
1723 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1724 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1725 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1726 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1727 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1728 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1729 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1730 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1731 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1732 cfmc write(iout,*)'RCUTINI=',rcutini
1733 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1734 cfmc write(iout,*)'GRAT=',grat
1735 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1736 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1740 c----------------------------------------------------------------------------
1742 implicit real*8 (a-h,o-z)
1743 include 'DIMENSIONS'
1744 include 'COMMON.MCM'
1745 include 'COMMON.MCE'
1746 include 'COMMON.IOUNITS'
1748 character*320 mcmcard
1749 call card_concat(mcmcard)
1750 call readi(mcmcard,'MAXACC',maxacc,100)
1751 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1752 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1753 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1754 call readi(mcmcard,'MAXREPM',maxrepm,200)
1755 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1756 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1757 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1758 call reada(mcmcard,'E_UP',e_up,5.0D0)
1759 call reada(mcmcard,'DELTE',delte,0.1D0)
1760 call readi(mcmcard,'NSWEEP',nsweep,5)
1761 call readi(mcmcard,'NSTEPH',nsteph,0)
1762 call readi(mcmcard,'NSTEPC',nstepc,0)
1763 call reada(mcmcard,'TMIN',tmin,298.0D0)
1764 call reada(mcmcard,'TMAX',tmax,298.0D0)
1765 call readi(mcmcard,'NWINDOW',nwindow,0)
1766 call readi(mcmcard,'PRINT_MC',print_mc,0)
1767 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1768 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1769 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1770 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1771 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1772 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1773 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1774 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1775 if (nwindow.gt.0) then
1776 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1778 winlen(i)=winend(i)-winstart(i)+1
1781 if (tmax.lt.tmin) tmax=tmin
1782 if (tmax.eq.tmin) then
1786 if (nstepc.gt.0 .and. nsteph.gt.0) then
1787 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1788 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1790 C Probabilities of different move types
1791 sumpro_type(0)=0.0D0
1792 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1793 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1794 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1795 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1796 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1797 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1798 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1800 print *,'i',i,' sumprotype',sumpro_type(i)
1801 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1802 print *,'i',i,' sumprotype',sumpro_type(i)
1806 c----------------------------------------------------------------------------
1807 subroutine read_minim
1808 implicit real*8 (a-h,o-z)
1809 include 'DIMENSIONS'
1810 include 'COMMON.MINIM'
1811 include 'COMMON.IOUNITS'
1813 character*320 minimcard
1814 call card_concat(minimcard)
1815 call readi(minimcard,'MAXMIN',maxmin,2000)
1816 call readi(minimcard,'MAXFUN',maxfun,5000)
1817 call readi(minimcard,'MINMIN',minmin,maxmin)
1818 call readi(minimcard,'MINFUN',minfun,maxmin)
1819 call reada(minimcard,'TOLF',tolf,1.0D-2)
1820 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1821 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1822 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1823 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1824 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1825 & 'Options in energy minimization:'
1826 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1827 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1828 & 'MinMin:',MinMin,' MinFun:',MinFun,
1829 & ' TolF:',TolF,' RTolF:',RTolF
1832 c----------------------------------------------------------------------------
1833 subroutine read_angles(kanal,*)
1834 implicit real*8 (a-h,o-z)
1835 include 'DIMENSIONS'
1836 include 'COMMON.GEO'
1837 include 'COMMON.VAR'
1838 include 'COMMON.CHAIN'
1839 include 'COMMON.IOUNITS'
1840 include 'COMMON.CONTROL'
1841 c Read angles from input
1843 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1844 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1845 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1846 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1849 c 9/7/01 avoid 180 deg valence angle
1850 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1852 theta(i)=deg2rad*theta(i)
1853 phi(i)=deg2rad*phi(i)
1854 alph(i)=deg2rad*alph(i)
1855 omeg(i)=deg2rad*omeg(i)
1860 c----------------------------------------------------------------------------
1861 subroutine reada(rekord,lancuch,wartosc,default)
1863 character*(*) rekord,lancuch
1864 double precision wartosc,default
1867 iread=index(rekord,lancuch)
1868 if (iread.eq.0) then
1872 iread=iread+ilen(lancuch)+1
1873 read (rekord(iread:),*,err=10,end=10) wartosc
1878 c----------------------------------------------------------------------------
1879 subroutine readi(rekord,lancuch,wartosc,default)
1881 character*(*) rekord,lancuch
1882 integer wartosc,default
1885 iread=index(rekord,lancuch)
1886 if (iread.eq.0) then
1890 iread=iread+ilen(lancuch)+1
1891 read (rekord(iread:),*,err=10,end=10) wartosc
1896 c----------------------------------------------------------------------------
1897 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1900 integer tablica(dim),default
1901 character*(*) rekord,lancuch
1908 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1909 if (iread.eq.0) return
1910 iread=iread+ilen(lancuch)+1
1911 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1914 c----------------------------------------------------------------------------
1915 subroutine multreada(rekord,lancuch,tablica,dim,default)
1918 double precision tablica(dim),default
1919 character*(*) rekord,lancuch
1926 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1927 if (iread.eq.0) return
1928 iread=iread+ilen(lancuch)+1
1929 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1932 c----------------------------------------------------------------------------
1933 subroutine openunits
1934 implicit real*8 (a-h,o-z)
1935 include 'DIMENSIONS'
1938 character*16 form,nodename
1941 include 'COMMON.SETUP'
1942 include 'COMMON.IOUNITS'
1944 include 'COMMON.CONTROL'
1945 integer lenpre,lenpot,ilen,lentmp
1947 character*3 out1file_text,ucase
1950 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1951 call getenv_loc("PREFIX",prefix)
1953 call getenv_loc("POT",pot)
1954 call getenv_loc("DIRTMP",tmpdir)
1955 call getenv_loc("CURDIR",curdir)
1956 call getenv_loc("OUT1FILE",out1file_text)
1957 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1958 out1file_text=ucase(out1file_text)
1959 if (out1file_text(1:1).eq."Y") then
1962 out1file=fg_rank.gt.0
1967 if (lentmp.gt.0) then
1968 write (*,'(80(1h!))')
1969 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1970 write (*,'(80(1h!))')
1971 write (*,*)"All output files will be on node /tmp directory."
1973 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1974 if (me.eq.king) then
1975 write (*,*) "The master node is ",nodename
1976 else if (fg_rank.eq.0) then
1977 write (*,*) "I am the CG slave node ",nodename
1979 write (*,*) "I am the FG slave node ",nodename
1982 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1983 lenpre = lentmp+lenpre+1
1985 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1986 C Get the names and open the input files
1987 #if defined(WINIFL) || defined(WINPGI)
1988 open(1,file=pref_orig(:ilen(pref_orig))//
1989 & '.inp',status='old',readonly,shared)
1990 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1991 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1992 C Get parameter filenames and open the parameter files.
1993 call getenv_loc('BONDPAR',bondname)
1994 open (ibond,file=bondname,status='old',readonly,shared)
1995 call getenv_loc('THETPAR',thetname)
1996 open (ithep,file=thetname,status='old',readonly,shared)
1997 call getenv_loc('ROTPAR',rotname)
1998 open (irotam,file=rotname,status='old',readonly,shared)
1999 call getenv_loc('TORPAR',torname)
2000 open (itorp,file=torname,status='old',readonly,shared)
2001 call getenv_loc('TORDPAR',tordname)
2002 open (itordp,file=tordname,status='old',readonly,shared)
2003 call getenv_loc('FOURIER',fouriername)
2004 open (ifourier,file=fouriername,status='old',readonly,shared)
2005 call getenv_loc('ELEPAR',elename)
2006 open (ielep,file=elename,status='old',readonly,shared)
2007 call getenv_loc('SIDEPAR',sidename)
2008 open (isidep,file=sidename,status='old',readonly,shared)
2009 #elif (defined CRAY) || (defined AIX)
2010 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2012 c print *,"Processor",myrank," opened file 1"
2013 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2014 c print *,"Processor",myrank," opened file 9"
2015 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2016 C Get parameter filenames and open the parameter files.
2017 call getenv_loc('BONDPAR',bondname)
2018 open (ibond,file=bondname,status='old',action='read')
2019 c print *,"Processor",myrank," opened file IBOND"
2020 call getenv_loc('THETPAR',thetname)
2021 open (ithep,file=thetname,status='old',action='read')
2022 c print *,"Processor",myrank," opened file ITHEP"
2023 call getenv_loc('ROTPAR',rotname)
2024 open (irotam,file=rotname,status='old',action='read')
2025 c print *,"Processor",myrank," opened file IROTAM"
2026 call getenv_loc('TORPAR',torname)
2027 open (itorp,file=torname,status='old',action='read')
2028 c print *,"Processor",myrank," opened file ITORP"
2029 call getenv_loc('TORDPAR',tordname)
2030 open (itordp,file=tordname,status='old',action='read')
2031 c print *,"Processor",myrank," opened file ITORDP"
2032 call getenv_loc('SCCORPAR',sccorname)
2033 open (isccor,file=sccorname,status='old',action='read')
2034 c print *,"Processor",myrank," opened file ISCCOR"
2035 call getenv_loc('FOURIER',fouriername)
2036 open (ifourier,file=fouriername,status='old',action='read')
2037 c print *,"Processor",myrank," opened file IFOURIER"
2038 call getenv_loc('ELEPAR',elename)
2039 open (ielep,file=elename,status='old',action='read')
2040 c print *,"Processor",myrank," opened file IELEP"
2041 call getenv_loc('SIDEPAR',sidename)
2042 open (isidep,file=sidename,status='old',action='read')
2043 c print *,"Processor",myrank," opened file ISIDEP"
2044 c print *,"Processor",myrank," opened parameter files"
2046 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2047 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2048 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2049 C Get parameter filenames and open the parameter files.
2050 call getenv_loc('BONDPAR',bondname)
2051 open (ibond,file=bondname,status='old')
2052 call getenv_loc('THETPAR',thetname)
2053 open (ithep,file=thetname,status='old')
2054 call getenv_loc('ROTPAR',rotname)
2055 open (irotam,file=rotname,status='old')
2056 call getenv_loc('TORPAR',torname)
2057 open (itorp,file=torname,status='old')
2058 call getenv_loc('TORDPAR',tordname)
2059 open (itordp,file=tordname,status='old')
2060 call getenv_loc('SCCORPAR',sccorname)
2061 open (isccor,file=sccorname,status='old')
2062 call getenv_loc('FOURIER',fouriername)
2063 open (ifourier,file=fouriername,status='old')
2064 call getenv_loc('ELEPAR',elename)
2065 open (ielep,file=elename,status='old')
2066 call getenv_loc('SIDEPAR',sidename)
2067 open (isidep,file=sidename,status='old')
2069 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2071 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2072 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2073 C Get parameter filenames and open the parameter files.
2074 call getenv_loc('BONDPAR',bondname)
2075 open (ibond,file=bondname,status='old',readonly)
2076 call getenv_loc('THETPAR',thetname)
2077 open (ithep,file=thetname,status='old',readonly)
2078 call getenv_loc('ROTPAR',rotname)
2079 open (irotam,file=rotname,status='old',readonly)
2080 call getenv_loc('TORPAR',torname)
2081 open (itorp,file=torname,status='old',readonly)
2082 call getenv_loc('TORDPAR',tordname)
2083 open (itordp,file=tordname,status='old',readonly)
2084 call getenv_loc('SCCORPAR',sccorname)
2085 open (isccor,file=sccorname,status='old',readonly)
2087 call getenv_loc('THETPARPDB',thetname_pdb)
2088 print *,"thetname_pdb ",thetname_pdb
2089 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2090 print *,ithep_pdb," opened"
2092 call getenv_loc('FOURIER',fouriername)
2093 open (ifourier,file=fouriername,status='old',readonly)
2094 call getenv_loc('ELEPAR',elename)
2095 open (ielep,file=elename,status='old',readonly)
2096 call getenv_loc('SIDEPAR',sidename)
2097 open (isidep,file=sidename,status='old',readonly)
2098 call getenv_loc('LIPTRANPAR',liptranname)
2099 open (iliptranpar,file=liptranname,status='old',action='read')
2101 call getenv_loc('ROTPARPDB',rotname_pdb)
2102 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2107 C 8/9/01 In the newest version SCp interaction constants are read from a file
2108 C Use -DOLDSCP to use hard-coded constants instead.
2110 call getenv_loc('SCPPAR',scpname)
2111 #if defined(WINIFL) || defined(WINPGI)
2112 open (iscpp,file=scpname,status='old',readonly,shared)
2113 #elif (defined CRAY) || (defined AIX)
2114 open (iscpp,file=scpname,status='old',action='read')
2116 open (iscpp,file=scpname,status='old')
2118 open (iscpp,file=scpname,status='old',readonly)
2121 call getenv_loc('PATTERN',patname)
2122 #if defined(WINIFL) || defined(WINPGI)
2123 open (icbase,file=patname,status='old',readonly,shared)
2124 #elif (defined CRAY) || (defined AIX)
2125 open (icbase,file=patname,status='old',action='read')
2127 open (icbase,file=patname,status='old')
2129 open (icbase,file=patname,status='old',readonly)
2132 C Open output file only for CG processes
2133 c print *,"Processor",myrank," fg_rank",fg_rank
2134 if (fg_rank.eq.0) then
2136 if (nodes.eq.1) then
2139 npos = dlog10(dfloat(nodes-1))+1
2141 if (npos.lt.3) npos=3
2142 write (liczba,'(i1)') npos
2143 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2145 write (liczba,form) me
2146 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2147 & liczba(:ilen(liczba))
2148 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2150 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2152 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2153 & liczba(:ilen(liczba))//'.mol2'
2154 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2155 & liczba(:ilen(liczba))//'.stat'
2157 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2158 & //liczba(:ilen(liczba))//'.stat')
2159 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2162 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2163 & liczba(:ilen(liczba))//'.const'
2168 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2169 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2170 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2171 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2172 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2174 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2176 rest2name=prefix(:ilen(prefix))//'.rst'
2178 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2181 #if defined(AIX) || defined(PGI)
2182 if (me.eq.king .or. .not. out1file)
2183 & open(iout,file=outname,status='unknown')
2185 if (fg_rank.gt.0) then
2186 write (liczba,'(i3.3)') myrank/nfgtasks
2187 write (ll,'(bz,i3.3)') fg_rank
2188 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2193 open(igeom,file=intname,status='unknown',position='append')
2194 open(ipdb,file=pdbname,status='unknown')
2195 open(imol2,file=mol2name,status='unknown')
2196 open(istat,file=statname,status='unknown',position='append')
2198 c1out open(iout,file=outname,status='unknown')
2201 if (me.eq.king .or. .not.out1file)
2202 & open(iout,file=outname,status='unknown')
2204 if (fg_rank.gt.0) then
2205 write (liczba,'(i3.3)') myrank/nfgtasks
2206 write (ll,'(bz,i3.3)') fg_rank
2207 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2212 open(igeom,file=intname,status='unknown',access='append')
2213 open(ipdb,file=pdbname,status='unknown')
2214 open(imol2,file=mol2name,status='unknown')
2215 open(istat,file=statname,status='unknown',access='append')
2217 c1out open(iout,file=outname,status='unknown')
2220 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2221 csa_seed=prefix(:lenpre)//'.CSA.seed'
2222 csa_history=prefix(:lenpre)//'.CSA.history'
2223 csa_bank=prefix(:lenpre)//'.CSA.bank'
2224 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2225 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2226 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2227 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2228 csa_int=prefix(:lenpre)//'.int'
2229 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2230 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2231 csa_in=prefix(:lenpre)//'.CSA.in'
2232 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2235 write (iout,'(80(1h-))')
2236 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2237 write (iout,'(80(1h-))')
2238 write (iout,*) "Input file : ",
2239 & pref_orig(:ilen(pref_orig))//'.inp'
2240 write (iout,*) "Output file : ",
2241 & outname(:ilen(outname))
2243 write (iout,*) "Sidechain potential file : ",
2244 & sidename(:ilen(sidename))
2246 write (iout,*) "SCp potential file : ",
2247 & scpname(:ilen(scpname))
2249 write (iout,*) "Electrostatic potential file : ",
2250 & elename(:ilen(elename))
2251 write (iout,*) "Cumulant coefficient file : ",
2252 & fouriername(:ilen(fouriername))
2253 write (iout,*) "Torsional parameter file : ",
2254 & torname(:ilen(torname))
2255 write (iout,*) "Double torsional parameter file : ",
2256 & tordname(:ilen(tordname))
2257 write (iout,*) "SCCOR parameter file : ",
2258 & sccorname(:ilen(sccorname))
2259 write (iout,*) "Bond & inertia constant file : ",
2260 & bondname(:ilen(bondname))
2261 write (iout,*) "Bending parameter file : ",
2262 & thetname(:ilen(thetname))
2263 write (iout,*) "Rotamer parameter file : ",
2264 & rotname(:ilen(rotname))
2265 write (iout,*) "Thetpdb parameter file : ",
2266 & thetname_pdb(:ilen(thetname_pdb))
2267 write (iout,*) "Threading database : ",
2268 & patname(:ilen(patname))
2270 &write (iout,*)" DIRTMP : ",
2272 write (iout,'(80(1h-))')
2276 c----------------------------------------------------------------------------
2277 subroutine card_concat(card)
2278 implicit real*8 (a-h,o-z)
2279 include 'DIMENSIONS'
2280 include 'COMMON.IOUNITS'
2282 character*80 karta,ucase
2284 read (inp,'(a)') karta
2287 do while (karta(80:80).eq.'&')
2288 card=card(:ilen(card)+1)//karta(:79)
2289 read (inp,'(a)') karta
2292 card=card(:ilen(card)+1)//karta
2295 c----------------------------------------------------------------------------------
2297 implicit real*8 (a-h,o-z)
2298 include 'DIMENSIONS'
2299 include 'COMMON.CHAIN'
2300 include 'COMMON.IOUNITS'
2302 open(irest2,file=rest2name,status='unknown')
2303 read(irest2,*) totT,EK,potE,totE,t_bath
2306 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2309 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2312 read (irest2,*) iset
2317 c---------------------------------------------------------------------------------
2318 subroutine read_fragments
2319 implicit real*8 (a-h,o-z)
2320 include 'DIMENSIONS'
2324 include 'COMMON.SETUP'
2325 include 'COMMON.CHAIN'
2326 include 'COMMON.IOUNITS'
2328 include 'COMMON.CONTROL'
2329 read(inp,*) nset,nfrag,npair,nfrag_back
2330 if(me.eq.king.or..not.out1file)
2331 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2332 & " nfrag_back",nfrag_back
2334 read(inp,*) mset(iset)
2336 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2338 if(me.eq.king.or..not.out1file)
2339 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2340 & ifrag(2,i,iset), qinfrag(i,iset)
2343 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2345 if(me.eq.king.or..not.out1file)
2346 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2347 & ipair(2,i,iset), qinpair(i,iset)
2350 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2351 & wfrag_back(3,i,iset),
2352 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2353 if(me.eq.king.or..not.out1file)
2354 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2355 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2360 C---------------------------------------------------------------------------
2361 subroutine read_afminp
2362 implicit real*8 (a-h,o-z)
2363 include 'DIMENSIONS'
2367 include 'COMMON.SETUP'
2368 include 'COMMON.CONTROL'
2369 include 'COMMON.CHAIN'
2370 include 'COMMON.IOUNITS'
2371 include 'COMMON.SBRIDGE'
2372 character*320 afmcard
2374 call card_concat(afmcard)
2375 call readi(afmcard,"BEG",afmbeg,0)
2376 call readi(afmcard,"END",afmend,0)
2377 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2378 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2379 print *,'FORCE=' ,forceAFMconst
2380 CCCC NOW PROPERTIES FOR AFM
2383 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2385 distafminit=dsqrt(distafminit)
2386 print *,'initdist',distafminit
2390 c-------------------------------------------------------------------------------
2391 subroutine read_dist_constr
2392 implicit real*8 (a-h,o-z)
2393 include 'DIMENSIONS'
2397 include 'COMMON.SETUP'
2398 include 'COMMON.CONTROL'
2399 include 'COMMON.CHAIN'
2400 include 'COMMON.IOUNITS'
2401 include 'COMMON.SBRIDGE'
2402 integer ifrag_(2,100),ipair_(2,100)
2403 double precision wfrag_(100),wpair_(100)
2404 character*500 controlcard
2406 write (iout,*) "Calling read_dist_constr"
2407 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2409 call card_concat(controlcard)
2410 call readi(controlcard,"NFRAG",nfrag_,0)
2411 call readi(controlcard,"NPAIR",npair_,0)
2412 call readi(controlcard,"NDIST",ndist_,0)
2413 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2414 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2415 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2416 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2417 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2418 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2419 c write (iout,*) "IFRAG"
2421 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2423 c write (iout,*) "IPAIR"
2425 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2429 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2430 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2431 & ifrag_(2,i)=nstart_sup+nsup-1
2432 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2434 if (wfrag_(i).gt.0.0d0) then
2435 do j=ifrag_(1,i),ifrag_(2,i)-1
2436 do k=j+1,ifrag_(2,i)
2437 c write (iout,*) "j",j," k",k
2439 if (constr_dist.eq.1) then
2444 forcon(nhpb)=wfrag_(i)
2445 else if (constr_dist.eq.2) then
2446 if (ddjk.le.dist_cut) then
2451 forcon(nhpb)=wfrag_(i)
2458 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2461 if (.not.out1file .or. me.eq.king)
2462 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2463 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2465 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2466 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2473 if (wpair_(i).gt.0.0d0) then
2481 do j=ifrag_(1,ii),ifrag_(2,ii)
2482 do k=ifrag_(1,jj),ifrag_(2,jj)
2486 forcon(nhpb)=wpair_(i)
2487 dhpb(nhpb)=dist(j,k)
2489 if (.not.out1file .or. me.eq.king)
2490 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2491 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2493 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2494 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2502 if (constr_dist.eq.11) then
2503 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2504 & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
2505 fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2508 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2509 & ibecarb(i),forcon(nhpb+1)
2511 if (forcon(nhpb+1).gt.0.0d0) then
2513 if (ibecarb(i).gt.0) then
2514 ihpb(i)=ihpb(i)+nres
2515 jhpb(i)=jhpb(i)+nres
2517 if (dhpb(nhpb).eq.0.0d0)
2518 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2520 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2521 C if (forcon(nhpb+1).gt.0.0d0) then
2523 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2525 if (.not.out1file .or. me.eq.king)
2526 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2527 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2529 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2530 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2537 c-------------------------------------------------------------------------------
2539 subroutine flush(iu)
2544 subroutine flush(iu)
2549 c------------------------------------------------------------------------------
2550 subroutine copy_to_tmp(source)
2551 include "DIMENSIONS"
2552 include "COMMON.IOUNITS"
2553 character*(*) source
2554 character* 256 tmpfile
2558 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2559 inquire(file=tmpfile,exist=ex)
2561 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2562 & " to temporary directory..."
2563 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2564 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2568 c------------------------------------------------------------------------------
2569 subroutine move_from_tmp(source)
2570 include "DIMENSIONS"
2571 include "COMMON.IOUNITS"
2572 character*(*) source
2575 write (*,*) "Moving ",source(:ilen(source)),
2576 & " from temporary directory to working directory"
2577 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2578 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2581 c------------------------------------------------------------------------------
2582 subroutine random_init(seed)
2584 C Initialize random number generator
2586 implicit real*8 (a-h,o-z)
2587 include 'DIMENSIONS'
2590 logical OKRandom, prng_restart
2592 integer iseed_array(4)
2594 include 'COMMON.IOUNITS'
2595 include 'COMMON.TIME1'
2596 include 'COMMON.THREAD'
2597 include 'COMMON.SBRIDGE'
2598 include 'COMMON.CONTROL'
2599 include 'COMMON.MCM'
2600 include 'COMMON.MAP'
2601 include 'COMMON.HEADER'
2602 include 'COMMON.CSA'
2603 include 'COMMON.CHAIN'
2604 include 'COMMON.MUCA'
2606 include 'COMMON.FFIELD'
2607 include 'COMMON.SETUP'
2608 iseed=-dint(dabs(seed))
2609 if (iseed.eq.0) then
2610 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2611 & 'Random seed undefined. The program will stop.'
2612 write (*,'(/80(1h*)/20x,a/80(1h*))')
2613 & 'Random seed undefined. The program will stop.'
2615 call mpi_finalize(mpi_comm_world,ierr)
2617 stop 'Bad random seed.'
2620 if (fg_rank.eq.0) then
2624 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2625 OKRandom = prng_restart(me,iseed)
2628 tmp=65536.0d0**(4-i)
2629 iseed_array(i) = dint(seed/tmp)
2630 seed=seed-iseed_array(i)*tmp
2633 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2634 & (iseed_array(i),i=1,4)
2635 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2636 & (iseed_array(i),i=1,4)
2637 OKRandom = prng_restart(me,iseed_array)
2640 c r1 = prng_next(me)
2641 r1=ran_number(0.0D0,1.0D0)
2643 & write (iout,*) 'ran_num',r1
2644 if (r1.lt.0.0d0) OKRandom=.false.
2646 if (.not.OKRandom) then
2647 write (iout,*) 'PRNG IS NOT WORKING!!!'
2648 print *,'PRNG IS NOT WORKING!!!'
2651 call mpi_abort(mpi_comm_world,error_msg,ierr)
2654 write (iout,*) 'too many processors for parallel prng'
2655 write (*,*) 'too many processors for parallel prng'
2663 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)