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,2800.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,'TORMODE',tor_mode,0)
159 C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
160 write(iout,*) "torsional and valence angle mode",tor_mode
161 call readi(controlcard,'MAXGEN',maxgen,10000)
162 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
163 call readi(controlcard,"KDIAG",kdiag,0)
164 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
165 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
166 & write (iout,*) "RESCALE_MODE",rescale_mode
167 split_ene=index(controlcard,'SPLIT_ENE').gt.0
168 if (index(controlcard,'REGULAR').gt.0.0D0) then
169 call reada(controlcard,'WEIDIS',weidis,0.1D0)
173 if (index(controlcard,'CHECKGRAD').gt.0) then
175 if (index(controlcard,'CART').gt.0) then
177 elseif (index(controlcard,'CARINT').gt.0) then
182 call reada(controlcard,'DELTA',aincr,1.0d-7)
183 elseif (index(controlcard,'THREAD').gt.0) then
185 call readi(controlcard,'THREAD',nthread,0)
186 if (nthread.gt.0) then
187 call reada(controlcard,'WEIDIS',weidis,0.1D0)
190 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
191 stop 'Error termination in Read_Control.'
193 else if (index(controlcard,'MCMA').gt.0) then
195 else if (index(controlcard,'MCEE').gt.0) then
197 else if (index(controlcard,'MULTCONF').gt.0) then
199 else if (index(controlcard,'MAP').gt.0) then
201 call readi(controlcard,'MAP',nmap,0)
202 else if (index(controlcard,'CSA').gt.0) then
204 crc else if (index(controlcard,'ZSCORE').gt.0) then
206 crc ZSCORE is rm from UNRES, modecalc=9 is available
209 cfcm else if (index(controlcard,'MCMF').gt.0) then
211 else if (index(controlcard,'SOFTREG').gt.0) then
213 else if (index(controlcard,'CHECK_BOND').gt.0) then
215 else if (index(controlcard,'TEST').gt.0) then
217 else if (index(controlcard,'MD').gt.0) then
219 else if (index(controlcard,'RE ').gt.0) then
223 lmuca=index(controlcard,'MUCA').gt.0
224 call readi(controlcard,'MUCADYN',mucadyn,0)
225 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
226 if (lmuca .and. (me.eq.king .or. .not.out1file ))
228 write (iout,*) 'MUCADYN=',mucadyn
229 write (iout,*) 'MUCASMOOTH=',muca_smooth
232 iscode=index(controlcard,'ONE_LETTER')
233 indphi=index(controlcard,'PHI')
234 indback=index(controlcard,'BACK')
235 iranconf=index(controlcard,'RAND_CONF')
236 i2ndstr=index(controlcard,'USE_SEC_PRED')
237 gradout=index(controlcard,'GRADOUT').gt.0
238 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
239 C DISTCHAINMAX become obsolete for periodic boundry condition
240 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
241 C Reading the dimensions of box in x,y,z coordinates
242 call reada(controlcard,'BOXX',boxxsize,100.0d0)
243 call reada(controlcard,'BOXY',boxysize,100.0d0)
244 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
245 c Cutoff range for interactions
246 call reada(controlcard,"R_CUT",r_cut,15.0d0)
247 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
248 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
249 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
250 if (lipthick.gt.0.0d0) then
251 bordliptop=(boxzsize+lipthick)/2.0
252 bordlipbot=bordliptop-lipthick
254 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
255 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
256 buflipbot=bordlipbot+lipbufthick
257 bufliptop=bordliptop-lipbufthick
258 if ((lipbufthick*2.0d0).gt.lipthick)
259 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
261 write(iout,*) "bordliptop=",bordliptop
262 write(iout,*) "bordlipbot=",bordlipbot
263 write(iout,*) "bufliptop=",bufliptop
264 write(iout,*) "buflipbot=",buflipbot
265 write (iout,*) "SHIELD MODE",shield_mode
266 if (shield_mode.gt.0) then
268 C VSolvSphere the volume of solving sphere
270 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
271 C there will be no distinction between proline peptide group and normal peptide
272 C group in case of shielding parameters
273 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
274 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
275 write (iout,*) VSolvSphere,VSolvSphere_div
276 C long axis of side chain
278 long_r_sidechain(i)=vbldsc0(1,i)
279 short_r_sidechain(i)=sigma0(i)
283 if (me.eq.king .or. .not.out1file )
284 & write (iout,*) "DISTCHAINMAX",distchainmax
286 if(me.eq.king.or..not.out1file)
287 & write (iout,'(2a)') diagmeth(kdiag),
288 & ' routine used to diagonalize matrices.'
291 c--------------------------------------------------------------------------
292 subroutine read_REMDpar
296 implicit real*8 (a-h,o-z)
298 include 'COMMON.IOUNITS'
299 include 'COMMON.TIME1'
302 include 'COMMON.LANGEVIN'
304 include 'COMMON.LANGEVIN.lang0'
306 include 'COMMON.INTERACT'
307 include 'COMMON.NAMES'
309 include 'COMMON.REMD'
310 include 'COMMON.CONTROL'
311 include 'COMMON.SETUP'
313 character*320 controlcard
314 character*3200 controlcard1
315 integer iremd_m_total
317 if(me.eq.king.or..not.out1file)
318 & write (iout,*) "REMD setup"
320 call card_concat(controlcard)
321 call readi(controlcard,"NREP",nrep,3)
322 call readi(controlcard,"NSTEX",nstex,1000)
323 call reada(controlcard,"RETMIN",retmin,10.0d0)
324 call reada(controlcard,"RETMAX",retmax,1000.0d0)
325 mremdsync=(index(controlcard,'SYNC').gt.0)
326 call readi(controlcard,"NSYN",i_sync_step,100)
327 restart1file=(index(controlcard,'REST1FILE').gt.0)
328 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
329 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
330 if(max_cache_traj_use.gt.max_cache_traj)
331 & max_cache_traj_use=max_cache_traj
332 if(me.eq.king.or..not.out1file) then
333 cd if (traj1file) then
334 crc caching is in testing - NTWX is not ignored
335 cd write (iout,*) "NTWX value is ignored"
336 cd write (iout,*) " trajectory is stored to one file by master"
337 cd write (iout,*) " before exchange at NSTEX intervals"
339 write (iout,*) "NREP= ",nrep
340 write (iout,*) "NSTEX= ",nstex
341 write (iout,*) "SYNC= ",mremdsync
342 write (iout,*) "NSYN= ",i_sync_step
343 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
346 if (index(controlcard,'TLIST').gt.0) then
348 call card_concat(controlcard1)
349 read(controlcard1,*) (remd_t(i),i=1,nrep)
350 if(me.eq.king.or..not.out1file)
351 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
354 if (index(controlcard,'MLIST').gt.0) then
356 call card_concat(controlcard1)
357 read(controlcard1,*) (remd_m(i),i=1,nrep)
358 if(me.eq.king.or..not.out1file) then
359 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
362 iremd_m_total=iremd_m_total+remd_m(i)
364 write (iout,*) 'Total number of replicas ',iremd_m_total
367 if(me.eq.king.or..not.out1file)
368 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
371 c--------------------------------------------------------------------------
372 subroutine read_MDpar
376 implicit real*8 (a-h,o-z)
378 include 'COMMON.IOUNITS'
379 include 'COMMON.TIME1'
382 include 'COMMON.LANGEVIN'
384 include 'COMMON.LANGEVIN.lang0'
386 include 'COMMON.INTERACT'
387 include 'COMMON.NAMES'
389 include 'COMMON.SETUP'
390 include 'COMMON.CONTROL'
391 include 'COMMON.SPLITELE'
393 character*320 controlcard
395 call card_concat(controlcard)
396 call readi(controlcard,"NSTEP",n_timestep,1000000)
397 call readi(controlcard,"NTWE",ntwe,100)
398 call readi(controlcard,"NTWX",ntwx,1000)
399 call reada(controlcard,"DT",d_time,1.0d-1)
400 call reada(controlcard,"DVMAX",dvmax,2.0d1)
401 call reada(controlcard,"DAMAX",damax,1.0d1)
402 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
403 call readi(controlcard,"LANG",lang,0)
404 RESPA = index(controlcard,"RESPA") .gt. 0
405 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
406 ntime_split0=ntime_split
407 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
408 ntime_split0=ntime_split
409 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
410 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
411 rest = index(controlcard,"REST").gt.0
412 tbf = index(controlcard,"TBF").gt.0
413 usampl = index(controlcard,"USAMPL").gt.0
415 mdpdb = index(controlcard,"MDPDB").gt.0
416 call reada(controlcard,"T_BATH",t_bath,300.0d0)
417 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
418 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
419 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
420 if (count_reset_moment.eq.0) count_reset_moment=1000000000
421 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
422 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
423 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
424 if (count_reset_vel.eq.0) count_reset_vel=1000000000
425 large = index(controlcard,"LARGE").gt.0
426 print_compon = index(controlcard,"PRINT_COMPON").gt.0
427 rattle = index(controlcard,"RATTLE").gt.0
428 c if performing umbrella sampling, fragments constrained are read from the fragment file
434 if(me.eq.king.or..not.out1file) then
436 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
438 write (iout,'(a)') "The units are:"
439 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
440 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
441 & " acceleration: angstrom/(48.9 fs)**2"
442 write (iout,'(a)') "energy: kcal/mol, temperature: K"
444 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
445 write (iout,'(a60,f10.5,a)')
446 & "Initial time step of numerical integration:",d_time,
448 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
450 write (iout,'(2a,i4,a)')
451 & "A-MTS algorithm used; initial time step for fast-varying",
452 & " short-range forces split into",ntime_split," steps."
453 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
454 & r_cut," lambda",rlamb
456 write (iout,'(2a,f10.5)')
457 & "Maximum acceleration threshold to reduce the time step",
458 & "/increase split number:",damax
459 write (iout,'(2a,f10.5)')
460 & "Maximum predicted energy drift to reduce the timestep",
461 & "/increase split number:",edriftmax
462 write (iout,'(a60,f10.5)')
463 & "Maximum velocity threshold to reduce velocities:",dvmax
464 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
465 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
466 if (rattle) write (iout,'(a60)')
467 & "Rattle algorithm used to constrain the virtual bonds"
471 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
472 call reada(controlcard,"RWAT",rwat,1.4d0)
473 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
474 surfarea=index(controlcard,"SURFAREA").gt.0
475 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
476 if(me.eq.king.or..not.out1file)then
477 write (iout,'(/a,$)') "Langevin dynamics calculation"
480 & " with direct integration of Langevin equations"
481 else if (lang.eq.2) then
482 write (iout,'(a/)') " with TINKER stochasic MD integrator"
483 else if (lang.eq.3) then
484 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
485 else if (lang.eq.4) then
486 write (iout,'(a/)') " in overdamped mode"
488 write (iout,'(//a,i5)')
489 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
492 write (iout,'(a60,f10.5)') "Temperature:",t_bath
493 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
494 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
495 write (iout,'(a60,f10.5)')
496 & "Scaling factor of the friction forces:",scal_fric
497 if (surfarea) write (iout,'(2a,i10,a)')
498 & "Friction coefficients will be scaled by solvent-accessible",
499 & " surface area every",reset_fricmat," steps."
501 c Calculate friction coefficients and bounds of stochastic forces
502 eta=6*pi*cPoise*etawat
503 if(me.eq.king.or..not.out1file)
504 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
506 gamp=scal_fric*(pstok+rwat)*eta
507 stdfp=dsqrt(2*Rb*t_bath/d_time)
509 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
510 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
512 if(me.eq.king.or..not.out1file)then
513 write (iout,'(/2a/)')
514 & "Radii of site types and friction coefficients and std's of",
515 & " stochastic forces of fully exposed sites"
516 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
518 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
519 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
523 if(me.eq.king.or..not.out1file)then
524 write (iout,'(a)') "Berendsen bath calculation"
525 write (iout,'(a60,f10.5)') "Temperature:",t_bath
526 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
528 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
529 & count_reset_moment," steps"
531 & write (iout,'(a,i10,a)')
532 & "Velocities will be reset at random every",count_reset_vel,
536 if(me.eq.king.or..not.out1file)
537 & write (iout,'(a31)') "Microcanonical mode calculation"
539 if(me.eq.king.or..not.out1file)then
540 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
542 write(iout,*) "MD running with constraints."
543 write(iout,*) "Equilibration time ", eq_time, " mtus."
544 write(iout,*) "Constraining ", nfrag," fragments."
545 write(iout,*) "Length of each fragment, weight and q0:"
547 write (iout,*) "Set of restraints #",iset
549 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
550 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
552 write(iout,*) "constraints between ", npair, "fragments."
553 write(iout,*) "constraint pairs, weights and q0:"
555 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
556 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
558 write(iout,*) "angle constraints within ", nfrag_back,
559 & "backbone fragments."
560 write(iout,*) "fragment, weights:"
562 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
563 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
564 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
567 iset=mod(kolor,nset)+1
570 if(me.eq.king.or..not.out1file)
571 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
574 c------------------------------------------------------------------------------
577 C Read molecular data.
579 implicit real*8 (a-h,o-z)
585 include 'COMMON.IOUNITS'
588 include 'COMMON.INTERACT'
589 include 'COMMON.LOCAL'
590 include 'COMMON.NAMES'
591 include 'COMMON.CHAIN'
592 include 'COMMON.FFIELD'
593 include 'COMMON.SBRIDGE'
594 include 'COMMON.HEADER'
595 include 'COMMON.CONTROL'
596 include 'COMMON.DBASE'
597 include 'COMMON.THREAD'
598 include 'COMMON.CONTACTS'
599 include 'COMMON.TORCNSTR'
600 include 'COMMON.TIME1'
601 include 'COMMON.BOUNDS'
603 include 'COMMON.SETUP'
604 include 'COMMON.SHIELD'
605 character*4 sequence(maxres)
607 double precision x(maxvar)
608 character*256 pdbfile
609 character*400 weightcard
610 character*80 weightcard_t,ucase
611 dimension itype_pdb(maxres)
612 common /pizda/ itype_pdb
613 logical seq_comp,fail
614 double precision energia(0:n_ene)
620 C Read weights of the subsequent energy terms.
621 call card_concat(weightcard)
622 call reada(weightcard,'WLONG',wlong,1.0D0)
623 call reada(weightcard,'WSC',wsc,wlong)
624 call reada(weightcard,'WSCP',wscp,wlong)
625 call reada(weightcard,'WELEC',welec,1.0D0)
626 call reada(weightcard,'WVDWPP',wvdwpp,welec)
627 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
628 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
629 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
630 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
631 call reada(weightcard,'WTURN3',wturn3,1.0D0)
632 call reada(weightcard,'WTURN4',wturn4,1.0D0)
633 call reada(weightcard,'WTURN6',wturn6,1.0D0)
634 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
635 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
636 call reada(weightcard,'WBOND',wbond,1.0D0)
637 call reada(weightcard,'WTOR',wtor,1.0D0)
638 call reada(weightcard,'WTORD',wtor_d,1.0D0)
639 call reada(weightcard,'WANG',wang,1.0D0)
640 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
641 call reada(weightcard,'SCAL14',scal14,0.4D0)
642 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
643 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
644 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
645 call reada(weightcard,'TEMP0',temp0,300.0d0)
646 call reada(weightcard,'WSHIELD',wshield,1.0d0)
647 call reada(weightcard,'WLT',wliptran,0.0D0)
648 if (index(weightcard,'SOFT').gt.0) ipot=6
649 C 12/1/95 Added weight for the multi-body term WCORR
650 call reada(weightcard,'WCORRH',wcorr,1.0D0)
651 if (wcorr4.gt.0.0d0) wcorr=wcorr4
671 if(me.eq.king.or..not.out1file)
672 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
673 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
675 10 format (/'Energy-term weights (unscaled):'//
676 & 'WSCC= ',f10.6,' (SC-SC)'/
677 & 'WSCP= ',f10.6,' (SC-p)'/
678 & 'WELEC= ',f10.6,' (p-p electr)'/
679 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
680 & 'WBOND= ',f10.6,' (stretching)'/
681 & 'WANG= ',f10.6,' (bending)'/
682 & 'WSCLOC= ',f10.6,' (SC local)'/
683 & 'WTOR= ',f10.6,' (torsional)'/
684 & 'WTORD= ',f10.6,' (double torsional)'/
685 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
686 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
687 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
688 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
689 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
690 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
691 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
692 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
693 & 'WTURN6= ',f10.6,' (turns, 6th order)')
694 if(me.eq.king.or..not.out1file)then
695 if (wcorr4.gt.0.0d0) then
696 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
697 & 'between contact pairs of peptide groups'
698 write (iout,'(2(a,f5.3/))')
699 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
700 & 'Range of quenching the correlation terms:',2*delt_corr
701 else if (wcorr.gt.0.0d0) then
702 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
703 & 'between contact pairs of peptide groups'
705 write (iout,'(a,f8.3)')
706 & 'Scaling factor of 1,4 SC-p interactions:',scal14
707 write (iout,'(a,f8.3)')
708 & 'General scaling factor of SC-p interactions:',scalscp
710 r0_corr=cutoff_corr-delt_corr
712 aad(i,1)=scalscp*aad(i,1)
713 aad(i,2)=scalscp*aad(i,2)
714 bad(i,1)=scalscp*bad(i,1)
715 bad(i,2)=scalscp*bad(i,2)
717 call rescale_weights(t_bath)
718 if(me.eq.king.or..not.out1file)
719 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
720 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
722 22 format (/'Energy-term weights (scaled):'//
723 & 'WSCC= ',f10.6,' (SC-SC)'/
724 & 'WSCP= ',f10.6,' (SC-p)'/
725 & 'WELEC= ',f10.6,' (p-p electr)'/
726 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
727 & 'WBOND= ',f10.6,' (stretching)'/
728 & 'WANG= ',f10.6,' (bending)'/
729 & 'WSCLOC= ',f10.6,' (SC local)'/
730 & 'WTOR= ',f10.6,' (torsional)'/
731 & 'WTORD= ',f10.6,' (double torsional)'/
732 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
733 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
734 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
735 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
736 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
737 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
738 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
739 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
740 & 'WTURN6= ',f10.6,' (turns, 6th order)')
741 if(me.eq.king.or..not.out1file)
742 & write (iout,*) "Reference temperature for weights calculation:",
744 call reada(weightcard,"D0CM",d0cm,3.78d0)
745 call reada(weightcard,"AKCM",akcm,15.1d0)
746 call reada(weightcard,"AKTH",akth,11.0d0)
747 call reada(weightcard,"AKCT",akct,12.0d0)
748 call reada(weightcard,"V1SS",v1ss,-1.08d0)
749 call reada(weightcard,"V2SS",v2ss,7.61d0)
750 call reada(weightcard,"V3SS",v3ss,13.7d0)
751 call reada(weightcard,"EBR",ebr,-5.50D0)
752 call reada(weightcard,"ATRISS",atriss,0.301D0)
753 call reada(weightcard,"BTRISS",btriss,0.021D0)
754 call reada(weightcard,"CTRISS",ctriss,1.001D0)
755 call reada(weightcard,"DTRISS",dtriss,1.001D0)
756 write (iout,*) "ATRISS=", atriss
757 write (iout,*) "BTRISS=", btriss
758 write (iout,*) "CTRISS=", ctriss
759 write (iout,*) "DTRISS=", dtriss
760 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
762 dyn_ss_mask(i)=.false.
766 dyn_ssbond_ij(i,j)=1.0d300
769 call reada(weightcard,"HT",Ht,0.0D0)
771 ss_depth=ebr/wsc-0.25*eps(1,1)
772 Ht=Ht/wsc-0.25*eps(1,1)
773 akcm=akcm*wstrain/wsc
774 akth=akth*wstrain/wsc
775 akct=akct*wstrain/wsc
776 v1ss=v1ss*wstrain/wsc
777 v2ss=v2ss*wstrain/wsc
778 v3ss=v3ss*wstrain/wsc
780 if (wstrain.ne.0.0) then
781 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
786 write (iout,*) "wshield,", wshield
787 if(me.eq.king.or..not.out1file) then
788 write (iout,*) "Parameters of the SS-bond potential:"
789 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
791 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
792 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
793 write (iout,*)" HT",Ht
794 print *,'indpdb=',indpdb,' pdbref=',pdbref
796 if (indpdb.gt.0 .or. pdbref) then
797 read(inp,'(a)') pdbfile
798 if(me.eq.king.or..not.out1file)
799 & write (iout,'(2a)') 'PDB data will be read from file ',
800 & pdbfile(:ilen(pdbfile))
801 open(ipdbin,file=pdbfile,status='old',err=33)
803 33 write (iout,'(a)') 'Error opening PDB file.'
806 c write (iout,*) 'Begin reading pdb data'
809 c write (iout,*) 'Finished reading pdb data'
811 if(me.eq.king.or..not.out1file)
812 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
813 & ' nstart_sup=',nstart_sup
815 itype_pdb(i)=itype(i)
819 nct=nstart_sup+nsup-1
820 call contact(.false.,ncont_ref,icont_ref,co)
823 if(me.eq.king.or..not.out1file)
824 & write(iout,*)'Adding sidechains'
828 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
831 do while (fail.and.nsi.le.maxsi)
832 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
835 if(fail) write(iout,*)'Adding sidechain failed for res ',
836 & i,' after ',nsi,' trials'
841 if (indpdb.eq.0) then
842 C Read sequence if not taken from the pdb file.
844 c print *,'nres=',nres
845 if (iscode.gt.0) then
846 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
848 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
850 C Convert sequence to numeric code
852 itype(i)=rescode(i,sequence(i),iscode)
854 C Assign initial virtual bond lengths
860 vbld(i+nres)=dsc(iabs(itype(i)))
861 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
862 c write (iout,*) "i",i," itype",itype(i),
863 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
867 c print '(20i4)',(itype(i),i=1,nres)
870 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
872 if (itype(i).eq.ntyp1) then
876 else if (iabs(itype(i+1)).ne.20) then
878 else if (iabs(itype(i)).ne.20) then
885 if(me.eq.king.or..not.out1file)then
886 write (iout,*) "ITEL"
888 write (iout,*) i,itype(i),itel(i)
890 print *,'Call Read_Bridge.'
893 C 8/13/98 Set limits to generating the dihedral angles
898 read (inp,*) ndih_constr
899 if (ndih_constr.gt.0) then
901 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
903 if(me.eq.king.or..not.out1file)then
905 & 'There are',ndih_constr,' constraints on phi angles.'
907 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
912 phi0(i)=deg2rad*phi0(i)
913 drange(i)=deg2rad*drange(i)
915 C if(me.eq.king.or..not.out1file)
916 C & write (iout,*) 'FTORS',ftors
919 phibound(1,ii) = phi0(i)-drange(i)
920 phibound(2,ii) = phi0(i)+drange(i)
923 C first setting the theta boundaries to 0 to pi
924 C this mean that there is no energy penalty for any angle occuring this can be applied
925 C for generate random conformation but is not implemented in this way
930 C begin reading theta constrains this is quartic constrains allowing to
931 C have smooth second derivative
932 if (with_theta_constr) then
933 C with_theta_constr is keyword allowing for occurance of theta constrains
934 read (inp,*) ntheta_constr
935 C ntheta_constr is the number of theta constrains
936 if (ntheta_constr.gt.0) then
938 read (inp,*) (itheta_constr(i),theta_constr0(i),
939 & theta_drange(i),for_thet_constr(i),
941 C the above code reads from 1 to ntheta_constr
942 C itheta_constr(i) residue i for which is theta_constr
943 C theta_constr0 the global minimum value
944 C theta_drange is range for which there is no energy penalty
945 C for_thet_constr is the force constant for quartic energy penalty
947 if(me.eq.king.or..not.out1file)then
949 & 'There are',ntheta_constr,' constraints on phi angles.'
951 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
957 theta_constr0(i)=deg2rad*theta_constr0(i)
958 theta_drange(i)=deg2rad*theta_drange(i)
960 C if(me.eq.king.or..not.out1file)
961 C & write (iout,*) 'FTORS',ftors
962 C do i=1,ntheta_constr
963 C ii = itheta_constr(i)
964 C thetabound(1,ii) = phi0(i)-drange(i)
965 C thetabound(2,ii) = phi0(i)+drange(i)
967 endif ! ntheta_constr.gt.0
968 endif! with_theta_constr
970 C with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
971 C write (iout,*) "with_dihed_constr ",with_dihed_constr
976 write (iout,'(a)') 'Boundaries in phi angle sampling:'
978 write (iout,'(a3,i5,2f10.1)')
979 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
985 cd print *,'NNT=',NNT,' NCT=',NCT
986 if (itype(1).eq.ntyp1) nnt=2
987 if (itype(nres).eq.ntyp1) nct=nct-1
989 if(me.eq.king.or..not.out1file)
990 & write (iout,'(a,i3)') 'nsup=',nsup
992 if (nsup.le.(nct-nnt+1)) then
993 do i=0,nct-nnt+1-nsup
994 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1000 & 'Error - sequences to be superposed do not match.'
1003 do i=0,nsup-(nct-nnt+1)
1004 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1006 nstart_sup=nstart_sup+i
1012 & 'Error - sequences to be superposed do not match.'
1015 if (nsup.eq.0) nsup=nct-nnt
1016 if (nstart_sup.eq.0) nstart_sup=nnt
1017 if (nstart_seq.eq.0) nstart_seq=nnt
1018 if(me.eq.king.or..not.out1file)
1019 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1020 & ' nstart_seq=',nstart_seq
1022 c--- Zscore rms -------
1023 if (nz_start.eq.0) nz_start=nnt
1024 if (nz_end.eq.0 .and. nsup.gt.0) then
1026 else if (nz_end.eq.0) then
1029 if(me.eq.king.or..not.out1file)then
1030 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1031 write (iout,*) 'IZ_SC=',iz_sc
1033 c----------------------
1036 if (.not.pdbref) then
1037 call read_angles(inp,*38)
1039 38 write (iout,'(a)') 'Error reading reference structure.'
1041 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1042 stop 'Error reading reference structure'
1044 39 call chainbuild_extconf
1046 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1053 cref(j,i,kkk)=c(j,i)
1056 call contact(.true.,ncont_ref,icont_ref,co)
1060 write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1062 if (constr_dist.gt.0) call read_dist_constr
1063 write (iout,*) "After read_dist_constr nhpb",nhpb
1064 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1066 if(me.eq.king.or..not.out1file)
1067 & write (iout,*) 'Contact order:',co
1069 if(me.eq.king.or..not.out1file)
1070 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1073 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1075 if(me.eq.king.or..not.out1file)
1076 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1077 & icont_ref(1,i),' ',
1078 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1082 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1083 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1084 & modecalc.ne.10) then
1085 C If input structure hasn't been supplied from the PDB file read or generate
1087 if (iranconf.eq.0 .and. .not. extconf) then
1088 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1089 & write (iout,'(a)') 'Initial geometry will be read in.'
1091 read(inp,'(8f10.5)',end=36,err=36)
1092 & ((c(l,k),l=1,3),k=1,nres),
1093 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1094 write (iout,*) "Exit READ_CART"
1095 write (iout,'(8f10.5)')
1096 & ((c(l,k),l=1,3),k=1,nres),
1097 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1098 call int_from_cart1(.true.)
1099 write (iout,*) "Finish INT_TO_CART"
1102 dc(j,i)=c(j,i+1)-c(j,i)
1103 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1107 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1109 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1110 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1116 call read_angles(inp,*36)
1117 call chainbuild_extconf
1120 36 write (iout,'(a)') 'Error reading angle file.'
1122 call mpi_finalize( MPI_COMM_WORLD,IERR )
1124 stop 'Error reading angle file.'
1126 else if (extconf) then
1127 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1128 & write (iout,'(a)') 'Extended chain initial geometry.'
1130 theta(i)=90d0*deg2rad
1133 phi(i)=180d0*deg2rad
1136 alph(i)=110d0*deg2rad
1139 omeg(i)=-120d0*deg2rad
1140 if (itype(i).le.0) omeg(i)=-omeg(i)
1142 call chainbuild_extconf
1144 if(me.eq.king.or..not.out1file)
1145 & write (iout,'(a)') 'Random-generated initial geometry.'
1149 if (me.eq.king .or. fg_rank.eq.0 .and. (
1150 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1154 call gen_rand_conf(itmp,*30)
1156 30 write (iout,*) 'Failed to generate random conformation',
1157 & ', itrial=',itrial
1158 write (*,*) 'Processor:',me,
1159 & ' Failed to generate random conformation',
1169 write (iout,'(a,i3,a)') 'Processor:',me,
1170 & ' error in generating random conformation.'
1171 write (*,'(a,i3,a)') 'Processor:',me,
1172 & ' error in generating random conformation.'
1175 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1180 & ' error in generating random conformation.'
1185 elseif (modecalc.eq.4) then
1186 read (inp,'(a)') intinname
1187 open (intin,file=intinname,status='old',err=333)
1188 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1189 & write (iout,'(a)') 'intinname',intinname
1190 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1192 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1194 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1196 stop 'Error opening angle file.'
1200 C Generate distance constraints, if the PDB structure is to be regularized.
1201 if (nthread.gt.0) then
1202 call read_threadbase
1205 if (me.eq.king .or. .not. out1file)
1207 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1208 write (iout,'(/a,i3,a)')
1209 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1210 write (iout,'(20i4)') (iss(i),i=1,ns)
1212 write(iout,*)"Running with dynamic disulfide-bond formation"
1214 write (iout,'(/a/)') 'Pre-formed links are:'
1220 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1221 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1227 if (ns.gt.0.and.dyn_ss) then
1231 forcon(i-nss)=forcon(i)
1238 dyn_ss_mask(iss(i))=.true.
1241 if (i2ndstr.gt.0) call secstrp2dihc
1242 c call geom_to_var(nvar,x)
1243 c call etotal(energia(0))
1244 c call enerprint(energia(0))
1245 c call briefout(0,etot)
1247 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1248 cd write (iout,'(a)') 'Variable list:'
1249 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1251 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1252 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1253 & 'Processor',myrank,': end reading molecular data.'
1258 c--------------------------------------------------------------------------
1259 logical function seq_comp(itypea,itypeb,length)
1261 integer length,itypea(length),itypeb(length)
1264 if (itypea(i).ne.itypeb(i)) then
1272 c-----------------------------------------------------------------------------
1273 subroutine read_bridge
1274 C Read information about disulfide bridges.
1275 implicit real*8 (a-h,o-z)
1276 include 'DIMENSIONS'
1280 include 'COMMON.IOUNITS'
1281 include 'COMMON.GEO'
1282 include 'COMMON.VAR'
1283 include 'COMMON.INTERACT'
1284 include 'COMMON.LOCAL'
1285 include 'COMMON.NAMES'
1286 include 'COMMON.CHAIN'
1287 include 'COMMON.FFIELD'
1288 include 'COMMON.SBRIDGE'
1289 include 'COMMON.HEADER'
1290 include 'COMMON.CONTROL'
1291 include 'COMMON.DBASE'
1292 include 'COMMON.THREAD'
1293 include 'COMMON.TIME1'
1294 include 'COMMON.SETUP'
1295 C Read bridging residues.
1296 read (inp,*) ns,(iss(i),i=1,ns)
1298 if(me.eq.king.or..not.out1file)
1299 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1300 C Check whether the specified bridging residues are cystines.
1302 if (itype(iss(i)).ne.1) then
1303 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1304 & 'Do you REALLY think that the residue ',
1305 & restyp(itype(iss(i))),i,
1306 & ' can form a disulfide bridge?!!!'
1307 write (*,'(2a,i3,a)')
1308 & 'Do you REALLY think that the residue ',
1309 & restyp(itype(iss(i))),i,
1310 & ' can form a disulfide bridge?!!!'
1312 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1317 C Read preformed bridges.
1319 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1321 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1324 C Check if the residues involved in bridges are in the specified list of
1325 C bridging residues.
1328 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1329 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1330 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1331 & ' contains residues present in other pairs.'
1332 write (*,'(a,i3,a)') 'Disulfide pair',i,
1333 & ' contains residues present in other pairs.'
1335 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1341 if (ihpb(i).eq.iss(j)) goto 10
1343 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1346 if (jhpb(i).eq.iss(j)) goto 20
1348 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1354 ihpb(i)=ihpb(i)+nres
1355 jhpb(i)=jhpb(i)+nres
1361 c----------------------------------------------------------------------------
1362 subroutine read_x(kanal,*)
1363 implicit real*8 (a-h,o-z)
1364 include 'DIMENSIONS'
1365 include 'COMMON.GEO'
1366 include 'COMMON.VAR'
1367 include 'COMMON.CHAIN'
1368 include 'COMMON.IOUNITS'
1369 include 'COMMON.CONTROL'
1370 include 'COMMON.LOCAL'
1371 include 'COMMON.INTERACT'
1372 c Read coordinates from input
1374 read(kanal,'(8f10.5)',end=10,err=10)
1375 & ((c(l,k),l=1,3),k=1,nres),
1376 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1379 c(j,2*nres)=c(j,nres)
1381 call int_from_cart1(.false.)
1384 dc(j,i)=c(j,i+1)-c(j,i)
1385 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1389 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1391 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1392 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1400 c----------------------------------------------------------------------------
1401 subroutine read_threadbase
1402 implicit real*8 (a-h,o-z)
1403 include 'DIMENSIONS'
1404 include 'COMMON.IOUNITS'
1405 include 'COMMON.GEO'
1406 include 'COMMON.VAR'
1407 include 'COMMON.INTERACT'
1408 include 'COMMON.LOCAL'
1409 include 'COMMON.NAMES'
1410 include 'COMMON.CHAIN'
1411 include 'COMMON.FFIELD'
1412 include 'COMMON.SBRIDGE'
1413 include 'COMMON.HEADER'
1414 include 'COMMON.CONTROL'
1415 include 'COMMON.DBASE'
1416 include 'COMMON.THREAD'
1417 include 'COMMON.TIME1'
1418 C Read pattern database for threading.
1419 read (icbase,*) nseq
1421 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1422 & nres_base(2,i),nres_base(3,i)
1423 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1425 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1426 c & nres_base(2,i),nres_base(3,i)
1427 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1431 if (weidis.eq.0.0D0) weidis=0.1D0
1440 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1441 write (iout,'(a,i5)') 'nexcl: ',nexcl
1442 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1445 c------------------------------------------------------------------------------
1446 subroutine setup_var
1447 implicit real*8 (a-h,o-z)
1448 include 'DIMENSIONS'
1449 include 'COMMON.IOUNITS'
1450 include 'COMMON.GEO'
1451 include 'COMMON.VAR'
1452 include 'COMMON.INTERACT'
1453 include 'COMMON.LOCAL'
1454 include 'COMMON.NAMES'
1455 include 'COMMON.CHAIN'
1456 include 'COMMON.FFIELD'
1457 include 'COMMON.SBRIDGE'
1458 include 'COMMON.HEADER'
1459 include 'COMMON.CONTROL'
1460 include 'COMMON.DBASE'
1461 include 'COMMON.THREAD'
1462 include 'COMMON.TIME1'
1463 C Set up variable list.
1469 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1471 ialph(i,1)=nvar+nside
1475 if (indphi.gt.0) then
1477 else if (indback.gt.0) then
1482 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1485 c----------------------------------------------------------------------------
1486 subroutine gen_dist_constr
1487 C Generate CA distance constraints.
1488 implicit real*8 (a-h,o-z)
1489 include 'DIMENSIONS'
1490 include 'COMMON.IOUNITS'
1491 include 'COMMON.GEO'
1492 include 'COMMON.VAR'
1493 include 'COMMON.INTERACT'
1494 include 'COMMON.LOCAL'
1495 include 'COMMON.NAMES'
1496 include 'COMMON.CHAIN'
1497 include 'COMMON.FFIELD'
1498 include 'COMMON.SBRIDGE'
1499 include 'COMMON.HEADER'
1500 include 'COMMON.CONTROL'
1501 include 'COMMON.DBASE'
1502 include 'COMMON.THREAD'
1503 include 'COMMON.TIME1'
1504 dimension itype_pdb(maxres)
1505 common /pizda/ itype_pdb
1507 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1508 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1509 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1511 do i=nstart_sup,nstart_sup+nsup-1
1512 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1513 cd & ' seq_pdb', restyp(itype_pdb(i))
1514 do j=i+2,nstart_sup+nsup-1
1516 ihpb(nhpb)=i+nstart_seq-nstart_sup
1517 jhpb(nhpb)=j+nstart_seq-nstart_sup
1519 dhpb(nhpb)=dist(i,j)
1522 cd write (iout,'(a)') 'Distance constraints:'
1527 cd if (ii.gt.nres) then
1532 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1533 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1534 cd & dhpb(i),forcon(i)
1538 c----------------------------------------------------------------------------
1540 implicit real*8 (a-h,o-z)
1541 include 'DIMENSIONS'
1542 include 'COMMON.MAP'
1543 include 'COMMON.IOUNITS'
1544 character*3 angid(4) /'THE','PHI','ALP','OME'/
1545 character*80 mapcard,ucase
1547 read (inp,'(a)') mapcard
1548 mapcard=ucase(mapcard)
1549 if (index(mapcard,'PHI').gt.0) then
1551 else if (index(mapcard,'THE').gt.0) then
1553 else if (index(mapcard,'ALP').gt.0) then
1555 else if (index(mapcard,'OME').gt.0) then
1558 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1559 stop 'Error - illegal variable spec in MAP card.'
1561 call readi (mapcard,'RES1',res1(imap),0)
1562 call readi (mapcard,'RES2',res2(imap),0)
1563 if (res1(imap).eq.0) then
1564 res1(imap)=res2(imap)
1565 else if (res2(imap).eq.0) then
1566 res2(imap)=res1(imap)
1568 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1570 & 'Error - illegal definition of variable group in MAP.'
1571 stop 'Error - illegal definition of variable group in MAP.'
1573 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1574 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1575 call readi(mapcard,'NSTEP',nstep(imap),0)
1576 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1578 & 'Illegal boundary and/or step size specification in MAP.'
1579 stop 'Illegal boundary and/or step size specification in MAP.'
1584 c----------------------------------------------------------------------------
1586 implicit real*8 (a-h,o-z)
1587 include 'DIMENSIONS'
1588 include 'COMMON.IOUNITS'
1589 include 'COMMON.GEO'
1590 include 'COMMON.CSA'
1591 include 'COMMON.BANK'
1592 include 'COMMON.CONTROL'
1594 character*620 mcmcard
1595 call card_concat(mcmcard)
1597 call readi(mcmcard,'NCONF',nconf,50)
1598 call readi(mcmcard,'NADD',nadd,0)
1599 call readi(mcmcard,'JSTART',jstart,1)
1600 call readi(mcmcard,'JEND',jend,1)
1601 call readi(mcmcard,'NSTMAX',nstmax,500000)
1602 call readi(mcmcard,'N0',n0,1)
1603 call readi(mcmcard,'N1',n1,6)
1604 call readi(mcmcard,'N2',n2,4)
1605 call readi(mcmcard,'N3',n3,0)
1606 call readi(mcmcard,'N4',n4,0)
1607 call readi(mcmcard,'N5',n5,0)
1608 call readi(mcmcard,'N6',n6,10)
1609 call readi(mcmcard,'N7',n7,0)
1610 call readi(mcmcard,'N8',n8,0)
1611 call readi(mcmcard,'N9',n9,0)
1612 call readi(mcmcard,'N14',n14,0)
1613 call readi(mcmcard,'N15',n15,0)
1614 call readi(mcmcard,'N16',n16,0)
1615 call readi(mcmcard,'N17',n17,0)
1616 call readi(mcmcard,'N18',n18,0)
1618 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1620 call readi(mcmcard,'NDIFF',ndiff,2)
1621 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1622 call readi(mcmcard,'IS1',is1,1)
1623 call readi(mcmcard,'IS2',is2,8)
1624 call readi(mcmcard,'NRAN0',nran0,4)
1625 call readi(mcmcard,'NRAN1',nran1,2)
1626 call readi(mcmcard,'IRR',irr,1)
1627 call readi(mcmcard,'NSEED',nseed,20)
1628 call readi(mcmcard,'NTOTAL',ntotal,10000)
1629 call reada(mcmcard,'CUT1',cut1,2.0d0)
1630 call reada(mcmcard,'CUT2',cut2,5.0d0)
1631 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1632 call readi(mcmcard,'ICMAX',icmax,3)
1633 call readi(mcmcard,'IRESTART',irestart,0)
1634 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1637 call reada(mcmcard,'DELE',dele,20.0d0)
1638 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1639 call readi(mcmcard,'IREF',iref,0)
1640 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1641 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1642 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1643 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1644 write (iout,*) "NCONF_IN",nconf_in
1647 c----------------------------------------------------------------------------
1648 cfmc subroutine mcmfread
1649 cfmc implicit real*8 (a-h,o-z)
1650 cfmc include 'DIMENSIONS'
1651 cfmc include 'COMMON.MCMF'
1652 cfmc include 'COMMON.IOUNITS'
1653 cfmc include 'COMMON.GEO'
1654 cfmc character*80 ucase
1655 cfmc character*620 mcmcard
1656 cfmc call card_concat(mcmcard)
1658 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1659 cfmc write(iout,*)'MAXRANT=',maxrant
1660 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1661 cfmc write(iout,*)'MAXFAM=',maxfam
1662 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1663 cfmc write(iout,*)'NNET1=',nnet1
1664 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1665 cfmc write(iout,*)'NNET2=',nnet2
1666 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1667 cfmc write(iout,*)'NNET3=',nnet3
1668 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1669 cfmc write(iout,*)'ILASTT=',ilastt
1670 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1671 cfmc write(iout,*)'MAXSTR=',maxstr
1672 cfmc maxstr_f=maxstr/maxfam
1673 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1674 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1675 cfmc write(iout,*)'NMCMF=',nmcmf
1676 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1677 cfmc write(iout,*)'IFOCUS=',ifocus
1678 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1679 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1680 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1681 cfmc write(iout,*)'INTPRT=',intprt
1682 cfmc call readi(mcmcard,'IPRT',iprt,100)
1683 cfmc write(iout,*)'IPRT=',iprt
1684 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1685 cfmc write(iout,*)'IMAXTR=',imaxtr
1686 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1687 cfmc write(iout,*)'MAXEVEN=',maxeven
1688 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1689 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1690 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1691 cfmc write(iout,*)'INIMIN=',inimin
1692 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1693 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1694 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1695 cfmc write(iout,*)'NTHREAD=',nthread
1696 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1697 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1698 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1699 cfmc write(iout,*)'MAXPERT=',maxpert
1700 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1701 cfmc write(iout,*)'IRMSD=',irmsd
1702 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1703 cfmc write(iout,*)'DENEMIN=',denemin
1704 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1705 cfmc write(iout,*)'RCUT1S=',rcut1s
1706 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1707 cfmc write(iout,*)'RCUT1E=',rcut1e
1708 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1709 cfmc write(iout,*)'RCUT2S=',rcut2s
1710 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1711 cfmc write(iout,*)'RCUT2E=',rcut2e
1712 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1713 cfmc write(iout,*)'DPERT1=',d_pert1
1714 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1715 cfmc write(iout,*)'DPERT1A=',d_pert1a
1716 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1717 cfmc write(iout,*)'DPERT2=',d_pert2
1718 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1719 cfmc write(iout,*)'DPERT2A=',d_pert2a
1720 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1721 cfmc write(iout,*)'DPERT2B=',d_pert2b
1722 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1723 cfmc write(iout,*)'DPERT2C=',d_pert2c
1724 cfmc d_pert1=deg2rad*d_pert1
1725 cfmc d_pert1a=deg2rad*d_pert1a
1726 cfmc d_pert2=deg2rad*d_pert2
1727 cfmc d_pert2a=deg2rad*d_pert2a
1728 cfmc d_pert2b=deg2rad*d_pert2b
1729 cfmc d_pert2c=deg2rad*d_pert2c
1730 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1731 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1732 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1733 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1734 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1735 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1736 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1737 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1738 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1739 cfmc write(iout,*)'RCUTINI=',rcutini
1740 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1741 cfmc write(iout,*)'GRAT=',grat
1742 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1743 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1747 c----------------------------------------------------------------------------
1749 implicit real*8 (a-h,o-z)
1750 include 'DIMENSIONS'
1751 include 'COMMON.MCM'
1752 include 'COMMON.MCE'
1753 include 'COMMON.IOUNITS'
1755 character*320 mcmcard
1756 call card_concat(mcmcard)
1757 call readi(mcmcard,'MAXACC',maxacc,100)
1758 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1759 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1760 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1761 call readi(mcmcard,'MAXREPM',maxrepm,200)
1762 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1763 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1764 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1765 call reada(mcmcard,'E_UP',e_up,5.0D0)
1766 call reada(mcmcard,'DELTE',delte,0.1D0)
1767 call readi(mcmcard,'NSWEEP',nsweep,5)
1768 call readi(mcmcard,'NSTEPH',nsteph,0)
1769 call readi(mcmcard,'NSTEPC',nstepc,0)
1770 call reada(mcmcard,'TMIN',tmin,298.0D0)
1771 call reada(mcmcard,'TMAX',tmax,298.0D0)
1772 call readi(mcmcard,'NWINDOW',nwindow,0)
1773 call readi(mcmcard,'PRINT_MC',print_mc,0)
1774 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1775 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1776 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1777 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1778 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1779 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1780 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1781 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1782 if (nwindow.gt.0) then
1783 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1785 winlen(i)=winend(i)-winstart(i)+1
1788 if (tmax.lt.tmin) tmax=tmin
1789 if (tmax.eq.tmin) then
1793 if (nstepc.gt.0 .and. nsteph.gt.0) then
1794 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1795 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1797 C Probabilities of different move types
1798 sumpro_type(0)=0.0D0
1799 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1800 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1801 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1802 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1803 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1804 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1805 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1807 print *,'i',i,' sumprotype',sumpro_type(i)
1808 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1809 print *,'i',i,' sumprotype',sumpro_type(i)
1813 c----------------------------------------------------------------------------
1814 subroutine read_minim
1815 implicit real*8 (a-h,o-z)
1816 include 'DIMENSIONS'
1817 include 'COMMON.MINIM'
1818 include 'COMMON.IOUNITS'
1820 character*320 minimcard
1821 call card_concat(minimcard)
1822 call readi(minimcard,'MAXMIN',maxmin,2000)
1823 call readi(minimcard,'MAXFUN',maxfun,5000)
1824 call readi(minimcard,'MINMIN',minmin,maxmin)
1825 call readi(minimcard,'MINFUN',minfun,maxmin)
1826 call reada(minimcard,'TOLF',tolf,1.0D-2)
1827 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1828 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1829 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1830 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1831 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1832 & 'Options in energy minimization:'
1833 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1834 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1835 & 'MinMin:',MinMin,' MinFun:',MinFun,
1836 & ' TolF:',TolF,' RTolF:',RTolF
1839 c----------------------------------------------------------------------------
1840 subroutine read_angles(kanal,*)
1841 implicit real*8 (a-h,o-z)
1842 include 'DIMENSIONS'
1843 include 'COMMON.GEO'
1844 include 'COMMON.VAR'
1845 include 'COMMON.CHAIN'
1846 include 'COMMON.IOUNITS'
1847 include 'COMMON.CONTROL'
1848 c Read angles from input
1850 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1851 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1852 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1853 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1856 c 9/7/01 avoid 180 deg valence angle
1857 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1859 theta(i)=deg2rad*theta(i)
1860 phi(i)=deg2rad*phi(i)
1861 alph(i)=deg2rad*alph(i)
1862 omeg(i)=deg2rad*omeg(i)
1867 c----------------------------------------------------------------------------
1868 subroutine reada(rekord,lancuch,wartosc,default)
1870 character*(*) rekord,lancuch
1871 double precision wartosc,default
1874 iread=index(rekord,lancuch)
1875 if (iread.eq.0) then
1879 iread=iread+ilen(lancuch)+1
1880 read (rekord(iread:),*,err=10,end=10) wartosc
1885 c----------------------------------------------------------------------------
1886 subroutine readi(rekord,lancuch,wartosc,default)
1888 character*(*) rekord,lancuch
1889 integer wartosc,default
1892 iread=index(rekord,lancuch)
1893 if (iread.eq.0) then
1897 iread=iread+ilen(lancuch)+1
1898 read (rekord(iread:),*,err=10,end=10) wartosc
1903 c----------------------------------------------------------------------------
1904 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1907 integer tablica(dim),default
1908 character*(*) rekord,lancuch
1915 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1916 if (iread.eq.0) return
1917 iread=iread+ilen(lancuch)+1
1918 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1921 c----------------------------------------------------------------------------
1922 subroutine multreada(rekord,lancuch,tablica,dim,default)
1925 double precision tablica(dim),default
1926 character*(*) rekord,lancuch
1933 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1934 if (iread.eq.0) return
1935 iread=iread+ilen(lancuch)+1
1936 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1939 c----------------------------------------------------------------------------
1940 subroutine openunits
1941 implicit real*8 (a-h,o-z)
1942 include 'DIMENSIONS'
1945 character*16 form,nodename
1948 include 'COMMON.SETUP'
1949 include 'COMMON.IOUNITS'
1951 include 'COMMON.CONTROL'
1952 integer lenpre,lenpot,ilen,lentmp
1954 character*3 out1file_text,ucase
1957 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1958 call getenv_loc("PREFIX",prefix)
1960 call getenv_loc("POT",pot)
1961 call getenv_loc("DIRTMP",tmpdir)
1962 call getenv_loc("CURDIR",curdir)
1963 call getenv_loc("OUT1FILE",out1file_text)
1964 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1965 out1file_text=ucase(out1file_text)
1966 if (out1file_text(1:1).eq."Y") then
1969 out1file=fg_rank.gt.0
1974 if (lentmp.gt.0) then
1975 write (*,'(80(1h!))')
1976 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1977 write (*,'(80(1h!))')
1978 write (*,*)"All output files will be on node /tmp directory."
1980 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1981 if (me.eq.king) then
1982 write (*,*) "The master node is ",nodename
1983 else if (fg_rank.eq.0) then
1984 write (*,*) "I am the CG slave node ",nodename
1986 write (*,*) "I am the FG slave node ",nodename
1989 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1990 lenpre = lentmp+lenpre+1
1992 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1993 C Get the names and open the input files
1994 #if defined(WINIFL) || defined(WINPGI)
1995 open(1,file=pref_orig(:ilen(pref_orig))//
1996 & '.inp',status='old',readonly,shared)
1997 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1998 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1999 C Get parameter filenames and open the parameter files.
2000 call getenv_loc('BONDPAR',bondname)
2001 open (ibond,file=bondname,status='old',readonly,shared)
2002 call getenv_loc('THETPAR',thetname)
2003 open (ithep,file=thetname,status='old',readonly,shared)
2004 call getenv_loc('ROTPAR',rotname)
2005 open (irotam,file=rotname,status='old',readonly,shared)
2006 call getenv_loc('TORPAR',torname)
2007 open (itorp,file=torname,status='old',readonly,shared)
2008 call getenv_loc('TORDPAR',tordname)
2009 open (itordp,file=tordname,status='old',readonly,shared)
2010 call getenv_loc('TORKCC',torkccname)
2011 open (itorkcc,file=torkccname,status='old',readonly,shared)
2012 call getenv_loc('THETKCC',thetkccname)
2013 open (ithetkcc,file=thetkccname,status='old',readonly,shared)
2014 call getenv_loc('FOURIER',fouriername)
2015 open (ifourier,file=fouriername,status='old',readonly,shared)
2016 call getenv_loc('ELEPAR',elename)
2017 open (ielep,file=elename,status='old',readonly,shared)
2018 call getenv_loc('SIDEPAR',sidename)
2019 open (isidep,file=sidename,status='old',readonly,shared)
2020 call getenv_loc('LIPTRANPAR',liptranname)
2021 open (iliptranpar,file=liptranname,status='old',readonly,shared)
2022 #elif (defined CRAY) || (defined AIX)
2023 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2025 c print *,"Processor",myrank," opened file 1"
2026 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2027 c print *,"Processor",myrank," opened file 9"
2028 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2029 C Get parameter filenames and open the parameter files.
2030 call getenv_loc('BONDPAR',bondname)
2031 open (ibond,file=bondname,status='old',action='read')
2032 c print *,"Processor",myrank," opened file IBOND"
2033 call getenv_loc('THETPAR',thetname)
2034 open (ithep,file=thetname,status='old',action='read')
2035 c print *,"Processor",myrank," opened file ITHEP"
2036 call getenv_loc('ROTPAR',rotname)
2037 open (irotam,file=rotname,status='old',action='read')
2038 c print *,"Processor",myrank," opened file IROTAM"
2039 call getenv_loc('TORPAR',torname)
2040 open (itorp,file=torname,status='old',action='read')
2041 c print *,"Processor",myrank," opened file ITORP"
2042 call getenv_loc('TORDPAR',tordname)
2043 open (itordp,file=tordname,status='old',action='read')
2044 call getenv_loc('TORKCC',torkccname)
2045 open (itorkcc,file=torkccname,status='old',action='read')
2046 call getenv_loc('THETKCC',thetkccname)
2047 open (ithetkcc,file=thetkccname,status='old',action='read')
2048 c print *,"Processor",myrank," opened file ITORDP"
2049 call getenv_loc('SCCORPAR',sccorname)
2050 open (isccor,file=sccorname,status='old',action='read')
2051 c print *,"Processor",myrank," opened file ISCCOR"
2052 call getenv_loc('FOURIER',fouriername)
2053 open (ifourier,file=fouriername,status='old',action='read')
2054 c print *,"Processor",myrank," opened file IFOURIER"
2055 call getenv_loc('ELEPAR',elename)
2056 open (ielep,file=elename,status='old',action='read')
2057 c print *,"Processor",myrank," opened file IELEP"
2058 call getenv_loc('SIDEPAR',sidename)
2059 open (isidep,file=sidename,status='old',action='read')
2060 call getenv_loc('LIPTRANPAR',liptranname)
2061 open (iliptranpar,file=liptranname,status='old',action='read')
2062 c print *,"Processor",myrank," opened file ISIDEP"
2063 c print *,"Processor",myrank," opened parameter files"
2065 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2066 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2067 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2068 C Get parameter filenames and open the parameter files.
2069 call getenv_loc('BONDPAR',bondname)
2070 open (ibond,file=bondname,status='old')
2071 call getenv_loc('THETPAR',thetname)
2072 open (ithep,file=thetname,status='old')
2073 call getenv_loc('ROTPAR',rotname)
2074 open (irotam,file=rotname,status='old')
2075 call getenv_loc('TORPAR',torname)
2076 open (itorp,file=torname,status='old')
2077 call getenv_loc('TORDPAR',tordname)
2078 open (itordp,file=tordname,status='old')
2079 call getenv_loc('TORKCC',torkccname)
2080 open (itorkcc,file=torkccname,status='old')
2081 call getenv_loc('THETKCC',thetkccname)
2082 open (ithetkcc,file=thetkccname,status='old')
2083 call getenv_loc('SCCORPAR',sccorname)
2084 open (isccor,file=sccorname,status='old')
2085 call getenv_loc('FOURIER',fouriername)
2086 open (ifourier,file=fouriername,status='old')
2087 call getenv_loc('ELEPAR',elename)
2088 open (ielep,file=elename,status='old')
2089 call getenv_loc('SIDEPAR',sidename)
2090 open (isidep,file=sidename,status='old')
2091 call getenv_loc('LIPTRANPAR',liptranname)
2092 open (iliptranpar,file=liptranname,status='old')
2094 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2096 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2097 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2098 C Get parameter filenames and open the parameter files.
2099 call getenv_loc('BONDPAR',bondname)
2100 open (ibond,file=bondname,status='old',readonly)
2101 call getenv_loc('THETPAR',thetname)
2102 open (ithep,file=thetname,status='old',readonly)
2103 call getenv_loc('ROTPAR',rotname)
2104 open (irotam,file=rotname,status='old',readonly)
2105 call getenv_loc('TORPAR',torname)
2106 open (itorp,file=torname,status='old',readonly)
2107 call getenv_loc('TORDPAR',tordname)
2108 open (itordp,file=tordname,status='old',readonly)
2109 call getenv_loc('TORKCC',torkccname)
2110 open (itorkcc,file=torkccname,status='old',readonly)
2111 call getenv_loc('THETKCC',thetkccname)
2112 open (ithetkcc,file=thetkccname,status='old',readonly)
2113 call getenv_loc('SCCORPAR',sccorname)
2114 open (isccor,file=sccorname,status='old',readonly)
2116 call getenv_loc('THETPARPDB',thetname_pdb)
2117 print *,"thetname_pdb ",thetname_pdb
2118 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2119 print *,ithep_pdb," opened"
2121 call getenv_loc('FOURIER',fouriername)
2122 open (ifourier,file=fouriername,status='old',readonly)
2123 call getenv_loc('ELEPAR',elename)
2124 open (ielep,file=elename,status='old',readonly)
2125 call getenv_loc('SIDEPAR',sidename)
2126 open (isidep,file=sidename,status='old',readonly)
2127 call getenv_loc('LIPTRANPAR',liptranname)
2128 open (iliptranpar,file=liptranname,status='old',action='read')
2130 call getenv_loc('ROTPARPDB',rotname_pdb)
2131 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2136 C 8/9/01 In the newest version SCp interaction constants are read from a file
2137 C Use -DOLDSCP to use hard-coded constants instead.
2139 call getenv_loc('SCPPAR',scpname)
2140 #if defined(WINIFL) || defined(WINPGI)
2141 open (iscpp,file=scpname,status='old',readonly,shared)
2142 #elif (defined CRAY) || (defined AIX)
2143 open (iscpp,file=scpname,status='old',action='read')
2145 open (iscpp,file=scpname,status='old')
2147 open (iscpp,file=scpname,status='old',readonly)
2150 call getenv_loc('PATTERN',patname)
2151 #if defined(WINIFL) || defined(WINPGI)
2152 open (icbase,file=patname,status='old',readonly,shared)
2153 #elif (defined CRAY) || (defined AIX)
2154 open (icbase,file=patname,status='old',action='read')
2156 open (icbase,file=patname,status='old')
2158 open (icbase,file=patname,status='old',readonly)
2161 C Open output file only for CG processes
2162 c print *,"Processor",myrank," fg_rank",fg_rank
2163 if (fg_rank.eq.0) then
2165 if (nodes.eq.1) then
2168 npos = dlog10(dfloat(nodes-1))+1
2170 if (npos.lt.3) npos=3
2171 write (liczba,'(i1)') npos
2172 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2174 write (liczba,form) me
2175 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2176 & liczba(:ilen(liczba))
2177 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2179 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2181 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2182 & liczba(:ilen(liczba))//'.mol2'
2183 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2184 & liczba(:ilen(liczba))//'.stat'
2186 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2187 & //liczba(:ilen(liczba))//'.stat')
2188 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2191 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2192 & liczba(:ilen(liczba))//'.const'
2197 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2198 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2199 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2200 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2201 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2203 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2205 rest2name=prefix(:ilen(prefix))//'.rst'
2207 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2210 #if defined(AIX) || defined(PGI)
2211 if (me.eq.king .or. .not. out1file)
2212 & open(iout,file=outname,status='unknown')
2214 if (fg_rank.gt.0) then
2215 write (liczba,'(i3.3)') myrank/nfgtasks
2216 write (ll,'(bz,i3.3)') fg_rank
2217 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2222 open(igeom,file=intname,status='unknown',position='append')
2223 open(ipdb,file=pdbname,status='unknown')
2224 open(imol2,file=mol2name,status='unknown')
2225 open(istat,file=statname,status='unknown',position='append')
2227 c1out open(iout,file=outname,status='unknown')
2230 if (me.eq.king .or. .not.out1file)
2231 & open(iout,file=outname,status='unknown')
2233 if (fg_rank.gt.0) then
2234 write (liczba,'(i3.3)') myrank/nfgtasks
2235 write (ll,'(bz,i3.3)') fg_rank
2236 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2241 open(igeom,file=intname,status='unknown',access='append')
2242 open(ipdb,file=pdbname,status='unknown')
2243 open(imol2,file=mol2name,status='unknown')
2244 open(istat,file=statname,status='unknown',access='append')
2246 c1out open(iout,file=outname,status='unknown')
2249 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2250 csa_seed=prefix(:lenpre)//'.CSA.seed'
2251 csa_history=prefix(:lenpre)//'.CSA.history'
2252 csa_bank=prefix(:lenpre)//'.CSA.bank'
2253 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2254 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2255 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2256 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2257 csa_int=prefix(:lenpre)//'.int'
2258 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2259 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2260 csa_in=prefix(:lenpre)//'.CSA.in'
2261 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2264 write (iout,'(80(1h-))')
2265 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2266 write (iout,'(80(1h-))')
2267 write (iout,*) "Input file : ",
2268 & pref_orig(:ilen(pref_orig))//'.inp'
2269 write (iout,*) "Output file : ",
2270 & outname(:ilen(outname))
2272 write (iout,*) "Sidechain potential file : ",
2273 & sidename(:ilen(sidename))
2275 write (iout,*) "SCp potential file : ",
2276 & scpname(:ilen(scpname))
2278 write (iout,*) "Electrostatic potential file : ",
2279 & elename(:ilen(elename))
2280 write (iout,*) "Cumulant coefficient file : ",
2281 & fouriername(:ilen(fouriername))
2282 write (iout,*) "Torsional parameter file : ",
2283 & torname(:ilen(torname))
2284 write (iout,*) "Double torsional parameter file : ",
2285 & tordname(:ilen(tordname))
2286 write (iout,*) "SCCOR parameter file : ",
2287 & sccorname(:ilen(sccorname))
2288 write (iout,*) "Bond & inertia constant file : ",
2289 & bondname(:ilen(bondname))
2290 write (iout,*) "Bending parameter file : ",
2291 & thetname(:ilen(thetname))
2292 write (iout,*) "Rotamer parameter file : ",
2293 & rotname(:ilen(rotname))
2294 write (iout,*) "Thetpdb parameter file : ",
2295 & thetname_pdb(:ilen(thetname_pdb))
2296 write (iout,*) "Threading database : ",
2297 & patname(:ilen(patname))
2299 &write (iout,*)" DIRTMP : ",
2301 write (iout,'(80(1h-))')
2305 c----------------------------------------------------------------------------
2306 subroutine card_concat(card)
2307 implicit real*8 (a-h,o-z)
2308 include 'DIMENSIONS'
2309 include 'COMMON.IOUNITS'
2311 character*80 karta,ucase
2313 read (inp,'(a)') karta
2316 do while (karta(80:80).eq.'&')
2317 card=card(:ilen(card)+1)//karta(:79)
2318 read (inp,'(a)') karta
2321 card=card(:ilen(card)+1)//karta
2324 c----------------------------------------------------------------------------------
2326 implicit real*8 (a-h,o-z)
2327 include 'DIMENSIONS'
2328 include 'COMMON.CHAIN'
2329 include 'COMMON.IOUNITS'
2331 open(irest2,file=rest2name,status='unknown')
2332 read(irest2,*) totT,EK,potE,totE,t_bath
2335 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2338 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2341 read (irest2,*) iset
2346 c---------------------------------------------------------------------------------
2347 subroutine read_fragments
2348 implicit real*8 (a-h,o-z)
2349 include 'DIMENSIONS'
2353 include 'COMMON.SETUP'
2354 include 'COMMON.CHAIN'
2355 include 'COMMON.IOUNITS'
2357 include 'COMMON.CONTROL'
2358 read(inp,*) nset,nfrag,npair,nfrag_back
2359 if(me.eq.king.or..not.out1file)
2360 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2361 & " nfrag_back",nfrag_back
2363 read(inp,*) mset(iset)
2365 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2367 if(me.eq.king.or..not.out1file)
2368 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2369 & ifrag(2,i,iset), qinfrag(i,iset)
2372 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2374 if(me.eq.king.or..not.out1file)
2375 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2376 & ipair(2,i,iset), qinpair(i,iset)
2379 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2380 & wfrag_back(3,i,iset),
2381 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2382 if(me.eq.king.or..not.out1file)
2383 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2384 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2389 C---------------------------------------------------------------------------
2390 subroutine read_afminp
2391 implicit real*8 (a-h,o-z)
2392 include 'DIMENSIONS'
2396 include 'COMMON.SETUP'
2397 include 'COMMON.CONTROL'
2398 include 'COMMON.CHAIN'
2399 include 'COMMON.IOUNITS'
2400 include 'COMMON.SBRIDGE'
2401 character*320 afmcard
2403 call card_concat(afmcard)
2404 call readi(afmcard,"BEG",afmbeg,0)
2405 call readi(afmcard,"END",afmend,0)
2406 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2407 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2408 print *,'FORCE=' ,forceAFMconst
2409 CCCC NOW PROPERTIES FOR AFM
2412 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2414 distafminit=dsqrt(distafminit)
2415 print *,'initdist',distafminit
2419 c-------------------------------------------------------------------------------
2420 subroutine read_dist_constr
2421 implicit real*8 (a-h,o-z)
2422 include 'DIMENSIONS'
2426 include 'COMMON.SETUP'
2427 include 'COMMON.CONTROL'
2428 include 'COMMON.CHAIN'
2429 include 'COMMON.IOUNITS'
2430 include 'COMMON.SBRIDGE'
2431 integer ifrag_(2,100),ipair_(2,100)
2432 double precision wfrag_(100),wpair_(100)
2433 character*500 controlcard
2435 write (iout,*) "Calling read_dist_constr"
2436 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2438 call card_concat(controlcard)
2439 call readi(controlcard,"NFRAG",nfrag_,0)
2440 call readi(controlcard,"NPAIR",npair_,0)
2441 call readi(controlcard,"NDIST",ndist_,0)
2442 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2443 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2444 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2445 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2446 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2447 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2448 c write (iout,*) "IFRAG"
2450 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2452 c write (iout,*) "IPAIR"
2454 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2458 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2459 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2460 & ifrag_(2,i)=nstart_sup+nsup-1
2461 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2463 if (wfrag_(i).gt.0.0d0) then
2464 do j=ifrag_(1,i),ifrag_(2,i)-1
2465 do k=j+1,ifrag_(2,i)
2466 c write (iout,*) "j",j," k",k
2468 if (constr_dist.eq.1) then
2473 forcon(nhpb)=wfrag_(i)
2474 else if (constr_dist.eq.2) then
2475 if (ddjk.le.dist_cut) then
2480 forcon(nhpb)=wfrag_(i)
2487 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2490 if (.not.out1file .or. me.eq.king)
2491 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2492 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2494 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2495 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2502 if (wpair_(i).gt.0.0d0) then
2510 do j=ifrag_(1,ii),ifrag_(2,ii)
2511 do k=ifrag_(1,jj),ifrag_(2,jj)
2515 forcon(nhpb)=wpair_(i)
2516 dhpb(nhpb)=dist(j,k)
2518 if (.not.out1file .or. me.eq.king)
2519 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2520 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2522 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2523 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2531 if (constr_dist.eq.11) then
2532 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2533 & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
2534 fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2537 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2538 & ibecarb(i),forcon(nhpb+1)
2540 if (forcon(nhpb+1).gt.0.0d0) then
2542 if (ibecarb(i).gt.0) then
2543 ihpb(i)=ihpb(i)+nres
2544 jhpb(i)=jhpb(i)+nres
2546 if (dhpb(nhpb).eq.0.0d0)
2547 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2549 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2550 C if (forcon(nhpb+1).gt.0.0d0) then
2552 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2554 if (.not.out1file .or. me.eq.king)
2555 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2556 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2558 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2559 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2566 c-------------------------------------------------------------------------------
2568 subroutine flush(iu)
2573 subroutine flush(iu)
2578 c------------------------------------------------------------------------------
2579 subroutine copy_to_tmp(source)
2580 include "DIMENSIONS"
2581 include "COMMON.IOUNITS"
2582 character*(*) source
2583 character* 256 tmpfile
2587 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2588 inquire(file=tmpfile,exist=ex)
2590 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2591 & " to temporary directory..."
2592 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2593 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2597 c------------------------------------------------------------------------------
2598 subroutine move_from_tmp(source)
2599 include "DIMENSIONS"
2600 include "COMMON.IOUNITS"
2601 character*(*) source
2604 write (*,*) "Moving ",source(:ilen(source)),
2605 & " from temporary directory to working directory"
2606 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2607 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2610 c------------------------------------------------------------------------------
2611 subroutine random_init(seed)
2613 C Initialize random number generator
2615 implicit real*8 (a-h,o-z)
2616 include 'DIMENSIONS'
2619 logical OKRandom, prng_restart
2621 integer iseed_array(4)
2623 include 'COMMON.IOUNITS'
2624 include 'COMMON.TIME1'
2625 include 'COMMON.THREAD'
2626 include 'COMMON.SBRIDGE'
2627 include 'COMMON.CONTROL'
2628 include 'COMMON.MCM'
2629 include 'COMMON.MAP'
2630 include 'COMMON.HEADER'
2631 include 'COMMON.CSA'
2632 include 'COMMON.CHAIN'
2633 include 'COMMON.MUCA'
2635 include 'COMMON.FFIELD'
2636 include 'COMMON.SETUP'
2637 iseed=-dint(dabs(seed))
2638 if (iseed.eq.0) then
2639 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2640 & 'Random seed undefined. The program will stop.'
2641 write (*,'(/80(1h*)/20x,a/80(1h*))')
2642 & 'Random seed undefined. The program will stop.'
2644 call mpi_finalize(mpi_comm_world,ierr)
2646 stop 'Bad random seed.'
2649 if (fg_rank.eq.0) then
2653 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2654 OKRandom = prng_restart(me,iseed)
2657 tmp=65536.0d0**(4-i)
2658 iseed_array(i) = dint(seed/tmp)
2659 seed=seed-iseed_array(i)*tmp
2662 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2663 & (iseed_array(i),i=1,4)
2664 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2665 & (iseed_array(i),i=1,4)
2666 OKRandom = prng_restart(me,iseed_array)
2669 c r1 = prng_next(me)
2670 r1=ran_number(0.0D0,1.0D0)
2672 & write (iout,*) 'ran_num',r1
2673 if (r1.lt.0.0d0) OKRandom=.false.
2675 if (.not.OKRandom) then
2676 write (iout,*) 'PRNG IS NOT WORKING!!!'
2677 print *,'PRNG IS NOT WORKING!!!'
2680 call mpi_abort(mpi_comm_world,error_msg,ierr)
2683 write (iout,*) 'too many processors for parallel prng'
2684 write (*,*) 'too many processors for parallel prng'
2692 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)