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,'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 elseif (index(controlcard,'THREAD').gt.0) then
184 call readi(controlcard,'THREAD',nthread,0)
185 if (nthread.gt.0) then
186 call reada(controlcard,'WEIDIS',weidis,0.1D0)
189 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
190 stop 'Error termination in Read_Control.'
192 else if (index(controlcard,'MCMA').gt.0) then
194 else if (index(controlcard,'MCEE').gt.0) then
196 else if (index(controlcard,'MULTCONF').gt.0) then
198 else if (index(controlcard,'MAP').gt.0) then
200 call readi(controlcard,'MAP',nmap,0)
201 else if (index(controlcard,'CSA').gt.0) then
203 crc else if (index(controlcard,'ZSCORE').gt.0) then
205 crc ZSCORE is rm from UNRES, modecalc=9 is available
208 cfcm else if (index(controlcard,'MCMF').gt.0) then
210 else if (index(controlcard,'SOFTREG').gt.0) then
212 else if (index(controlcard,'CHECK_BOND').gt.0) then
214 else if (index(controlcard,'TEST').gt.0) then
216 else if (index(controlcard,'MD').gt.0) then
218 else if (index(controlcard,'RE ').gt.0) then
222 lmuca=index(controlcard,'MUCA').gt.0
223 call readi(controlcard,'MUCADYN',mucadyn,0)
224 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
225 if (lmuca .and. (me.eq.king .or. .not.out1file ))
227 write (iout,*) 'MUCADYN=',mucadyn
228 write (iout,*) 'MUCASMOOTH=',muca_smooth
231 iscode=index(controlcard,'ONE_LETTER')
232 indphi=index(controlcard,'PHI')
233 indback=index(controlcard,'BACK')
234 iranconf=index(controlcard,'RAND_CONF')
235 i2ndstr=index(controlcard,'USE_SEC_PRED')
236 gradout=index(controlcard,'GRADOUT').gt.0
237 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
238 C DISTCHAINMAX become obsolete for periodic boundry condition
239 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
240 C Reading the dimensions of box in x,y,z coordinates
241 call reada(controlcard,'BOXX',boxxsize,100.0d0)
242 call reada(controlcard,'BOXY',boxysize,100.0d0)
243 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
244 c Cutoff range for interactions
245 call reada(controlcard,"R_CUT",r_cut,15.0d0)
246 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
247 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
248 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
249 if (lipthick.gt.0.0d0) then
250 bordliptop=(boxzsize+lipthick)/2.0
251 bordlipbot=bordliptop-lipthick
253 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
254 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
255 buflipbot=bordlipbot+lipbufthick
256 bufliptop=bordliptop-lipbufthick
257 if ((lipbufthick*2.0d0).gt.lipthick)
258 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
260 write(iout,*) "bordliptop=",bordliptop
261 write(iout,*) "bordlipbot=",bordlipbot
262 write(iout,*) "bufliptop=",bufliptop
263 write(iout,*) "buflipbot=",buflipbot
264 write (iout,*) "SHIELD MODE",shield_mode
265 if (shield_mode.gt.0) then
267 C VSolvSphere the volume of solving sphere
269 C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
270 C there will be no distinction between proline peptide group and normal peptide
271 C group in case of shielding parameters
272 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
273 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
274 write (iout,*) VSolvSphere,VSolvSphere_div
275 C long axis of side chain
277 long_r_sidechain(i)=vbldsc0(1,i)
278 short_r_sidechain(i)=sigma0(i)
282 if (me.eq.king .or. .not.out1file )
283 & write (iout,*) "DISTCHAINMAX",distchainmax
285 if(me.eq.king.or..not.out1file)
286 & write (iout,'(2a)') diagmeth(kdiag),
287 & ' routine used to diagonalize matrices.'
290 c--------------------------------------------------------------------------
291 subroutine read_REMDpar
295 implicit real*8 (a-h,o-z)
297 include 'COMMON.IOUNITS'
298 include 'COMMON.TIME1'
301 include 'COMMON.LANGEVIN'
303 include 'COMMON.LANGEVIN.lang0'
305 include 'COMMON.INTERACT'
306 include 'COMMON.NAMES'
308 include 'COMMON.REMD'
309 include 'COMMON.CONTROL'
310 include 'COMMON.SETUP'
312 character*320 controlcard
313 character*3200 controlcard1
314 integer iremd_m_total
316 if(me.eq.king.or..not.out1file)
317 & write (iout,*) "REMD setup"
319 call card_concat(controlcard)
320 call readi(controlcard,"NREP",nrep,3)
321 call readi(controlcard,"NSTEX",nstex,1000)
322 call reada(controlcard,"RETMIN",retmin,10.0d0)
323 call reada(controlcard,"RETMAX",retmax,1000.0d0)
324 mremdsync=(index(controlcard,'SYNC').gt.0)
325 call readi(controlcard,"NSYN",i_sync_step,100)
326 restart1file=(index(controlcard,'REST1FILE').gt.0)
327 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
328 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
329 if(max_cache_traj_use.gt.max_cache_traj)
330 & max_cache_traj_use=max_cache_traj
331 if(me.eq.king.or..not.out1file) then
332 cd if (traj1file) then
333 crc caching is in testing - NTWX is not ignored
334 cd write (iout,*) "NTWX value is ignored"
335 cd write (iout,*) " trajectory is stored to one file by master"
336 cd write (iout,*) " before exchange at NSTEX intervals"
338 write (iout,*) "NREP= ",nrep
339 write (iout,*) "NSTEX= ",nstex
340 write (iout,*) "SYNC= ",mremdsync
341 write (iout,*) "NSYN= ",i_sync_step
342 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
345 if (index(controlcard,'TLIST').gt.0) then
347 call card_concat(controlcard1)
348 read(controlcard1,*) (remd_t(i),i=1,nrep)
349 if(me.eq.king.or..not.out1file)
350 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
353 if (index(controlcard,'MLIST').gt.0) then
355 call card_concat(controlcard1)
356 read(controlcard1,*) (remd_m(i),i=1,nrep)
357 if(me.eq.king.or..not.out1file) then
358 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
361 iremd_m_total=iremd_m_total+remd_m(i)
363 write (iout,*) 'Total number of replicas ',iremd_m_total
366 if(me.eq.king.or..not.out1file)
367 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
370 c--------------------------------------------------------------------------
371 subroutine read_MDpar
375 implicit real*8 (a-h,o-z)
377 include 'COMMON.IOUNITS'
378 include 'COMMON.TIME1'
381 include 'COMMON.LANGEVIN'
383 include 'COMMON.LANGEVIN.lang0'
385 include 'COMMON.INTERACT'
386 include 'COMMON.NAMES'
388 include 'COMMON.SETUP'
389 include 'COMMON.CONTROL'
390 include 'COMMON.SPLITELE'
392 character*320 controlcard
394 call card_concat(controlcard)
395 call readi(controlcard,"NSTEP",n_timestep,1000000)
396 call readi(controlcard,"NTWE",ntwe,100)
397 call readi(controlcard,"NTWX",ntwx,1000)
398 call reada(controlcard,"DT",d_time,1.0d-1)
399 call reada(controlcard,"DVMAX",dvmax,2.0d1)
400 call reada(controlcard,"DAMAX",damax,1.0d1)
401 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
402 call readi(controlcard,"LANG",lang,0)
403 RESPA = index(controlcard,"RESPA") .gt. 0
404 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
405 ntime_split0=ntime_split
406 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
407 ntime_split0=ntime_split
408 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
409 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
410 rest = index(controlcard,"REST").gt.0
411 tbf = index(controlcard,"TBF").gt.0
412 usampl = index(controlcard,"USAMPL").gt.0
414 mdpdb = index(controlcard,"MDPDB").gt.0
415 call reada(controlcard,"T_BATH",t_bath,300.0d0)
416 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
417 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
418 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
419 if (count_reset_moment.eq.0) count_reset_moment=1000000000
420 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
421 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
422 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
423 if (count_reset_vel.eq.0) count_reset_vel=1000000000
424 large = index(controlcard,"LARGE").gt.0
425 print_compon = index(controlcard,"PRINT_COMPON").gt.0
426 rattle = index(controlcard,"RATTLE").gt.0
427 c if performing umbrella sampling, fragments constrained are read from the fragment file
433 if(me.eq.king.or..not.out1file) then
435 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
437 write (iout,'(a)') "The units are:"
438 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
439 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
440 & " acceleration: angstrom/(48.9 fs)**2"
441 write (iout,'(a)') "energy: kcal/mol, temperature: K"
443 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
444 write (iout,'(a60,f10.5,a)')
445 & "Initial time step of numerical integration:",d_time,
447 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
449 write (iout,'(2a,i4,a)')
450 & "A-MTS algorithm used; initial time step for fast-varying",
451 & " short-range forces split into",ntime_split," steps."
452 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
453 & r_cut," lambda",rlamb
455 write (iout,'(2a,f10.5)')
456 & "Maximum acceleration threshold to reduce the time step",
457 & "/increase split number:",damax
458 write (iout,'(2a,f10.5)')
459 & "Maximum predicted energy drift to reduce the timestep",
460 & "/increase split number:",edriftmax
461 write (iout,'(a60,f10.5)')
462 & "Maximum velocity threshold to reduce velocities:",dvmax
463 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
464 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
465 if (rattle) write (iout,'(a60)')
466 & "Rattle algorithm used to constrain the virtual bonds"
470 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
471 call reada(controlcard,"RWAT",rwat,1.4d0)
472 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
473 surfarea=index(controlcard,"SURFAREA").gt.0
474 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
475 if(me.eq.king.or..not.out1file)then
476 write (iout,'(/a,$)') "Langevin dynamics calculation"
479 & " with direct integration of Langevin equations"
480 else if (lang.eq.2) then
481 write (iout,'(a/)') " with TINKER stochasic MD integrator"
482 else if (lang.eq.3) then
483 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
484 else if (lang.eq.4) then
485 write (iout,'(a/)') " in overdamped mode"
487 write (iout,'(//a,i5)')
488 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
491 write (iout,'(a60,f10.5)') "Temperature:",t_bath
492 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
493 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
494 write (iout,'(a60,f10.5)')
495 & "Scaling factor of the friction forces:",scal_fric
496 if (surfarea) write (iout,'(2a,i10,a)')
497 & "Friction coefficients will be scaled by solvent-accessible",
498 & " surface area every",reset_fricmat," steps."
500 c Calculate friction coefficients and bounds of stochastic forces
501 eta=6*pi*cPoise*etawat
502 if(me.eq.king.or..not.out1file)
503 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
505 gamp=scal_fric*(pstok+rwat)*eta
506 stdfp=dsqrt(2*Rb*t_bath/d_time)
508 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
509 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
511 if(me.eq.king.or..not.out1file)then
512 write (iout,'(/2a/)')
513 & "Radii of site types and friction coefficients and std's of",
514 & " stochastic forces of fully exposed sites"
515 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
517 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
518 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
522 if(me.eq.king.or..not.out1file)then
523 write (iout,'(a)') "Berendsen bath calculation"
524 write (iout,'(a60,f10.5)') "Temperature:",t_bath
525 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
527 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
528 & count_reset_moment," steps"
530 & write (iout,'(a,i10,a)')
531 & "Velocities will be reset at random every",count_reset_vel,
535 if(me.eq.king.or..not.out1file)
536 & write (iout,'(a31)') "Microcanonical mode calculation"
538 if(me.eq.king.or..not.out1file)then
539 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
541 write(iout,*) "MD running with constraints."
542 write(iout,*) "Equilibration time ", eq_time, " mtus."
543 write(iout,*) "Constraining ", nfrag," fragments."
544 write(iout,*) "Length of each fragment, weight and q0:"
546 write (iout,*) "Set of restraints #",iset
548 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
549 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
551 write(iout,*) "constraints between ", npair, "fragments."
552 write(iout,*) "constraint pairs, weights and q0:"
554 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
555 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
557 write(iout,*) "angle constraints within ", nfrag_back,
558 & "backbone fragments."
559 write(iout,*) "fragment, weights:"
561 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
562 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
563 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
566 iset=mod(kolor,nset)+1
569 if(me.eq.king.or..not.out1file)
570 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
573 c------------------------------------------------------------------------------
576 C Read molecular data.
578 implicit real*8 (a-h,o-z)
584 include 'COMMON.IOUNITS'
587 include 'COMMON.INTERACT'
588 include 'COMMON.LOCAL'
589 include 'COMMON.NAMES'
590 include 'COMMON.CHAIN'
591 include 'COMMON.FFIELD'
592 include 'COMMON.SBRIDGE'
593 include 'COMMON.HEADER'
594 include 'COMMON.CONTROL'
595 include 'COMMON.DBASE'
596 include 'COMMON.THREAD'
597 include 'COMMON.CONTACTS'
598 include 'COMMON.TORCNSTR'
599 include 'COMMON.TIME1'
600 include 'COMMON.BOUNDS'
602 include 'COMMON.SETUP'
603 character*4 sequence(maxres)
605 double precision x(maxvar)
606 character*256 pdbfile
607 character*400 weightcard
608 character*80 weightcard_t,ucase
609 dimension itype_pdb(maxres)
610 common /pizda/ itype_pdb
611 logical seq_comp,fail
612 double precision energia(0:n_ene)
618 C Read weights of the subsequent energy terms.
619 call card_concat(weightcard)
620 call reada(weightcard,'WLONG',wlong,1.0D0)
621 call reada(weightcard,'WSC',wsc,wlong)
622 call reada(weightcard,'WSCP',wscp,wlong)
623 call reada(weightcard,'WELEC',welec,1.0D0)
624 call reada(weightcard,'WVDWPP',wvdwpp,welec)
625 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
626 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
627 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
628 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
629 call reada(weightcard,'WTURN3',wturn3,1.0D0)
630 call reada(weightcard,'WTURN4',wturn4,1.0D0)
631 call reada(weightcard,'WTURN6',wturn6,1.0D0)
632 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
633 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
634 call reada(weightcard,'WBOND',wbond,1.0D0)
635 call reada(weightcard,'WTOR',wtor,1.0D0)
636 call reada(weightcard,'WTORD',wtor_d,1.0D0)
637 call reada(weightcard,'WANG',wang,1.0D0)
638 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
639 call reada(weightcard,'SCAL14',scal14,0.4D0)
640 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
641 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
642 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
643 call reada(weightcard,'TEMP0',temp0,300.0d0)
644 call reada(weightcard,'WLT',wliptran,0.0D0)
645 if (index(weightcard,'SOFT').gt.0) ipot=6
646 C 12/1/95 Added weight for the multi-body term WCORR
647 call reada(weightcard,'WCORRH',wcorr,1.0D0)
648 if (wcorr4.gt.0.0d0) wcorr=wcorr4
668 if(me.eq.king.or..not.out1file)
669 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
670 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
672 10 format (/'Energy-term weights (unscaled):'//
673 & 'WSCC= ',f10.6,' (SC-SC)'/
674 & 'WSCP= ',f10.6,' (SC-p)'/
675 & 'WELEC= ',f10.6,' (p-p electr)'/
676 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
677 & 'WBOND= ',f10.6,' (stretching)'/
678 & 'WANG= ',f10.6,' (bending)'/
679 & 'WSCLOC= ',f10.6,' (SC local)'/
680 & 'WTOR= ',f10.6,' (torsional)'/
681 & 'WTORD= ',f10.6,' (double torsional)'/
682 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
683 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
684 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
685 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
686 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
687 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
688 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
689 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
690 & 'WTURN6= ',f10.6,' (turns, 6th order)')
691 if(me.eq.king.or..not.out1file)then
692 if (wcorr4.gt.0.0d0) then
693 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
694 & 'between contact pairs of peptide groups'
695 write (iout,'(2(a,f5.3/))')
696 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
697 & 'Range of quenching the correlation terms:',2*delt_corr
698 else if (wcorr.gt.0.0d0) then
699 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
700 & 'between contact pairs of peptide groups'
702 write (iout,'(a,f8.3)')
703 & 'Scaling factor of 1,4 SC-p interactions:',scal14
704 write (iout,'(a,f8.3)')
705 & 'General scaling factor of SC-p interactions:',scalscp
707 r0_corr=cutoff_corr-delt_corr
709 aad(i,1)=scalscp*aad(i,1)
710 aad(i,2)=scalscp*aad(i,2)
711 bad(i,1)=scalscp*bad(i,1)
712 bad(i,2)=scalscp*bad(i,2)
714 call rescale_weights(t_bath)
715 if(me.eq.king.or..not.out1file)
716 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
717 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
719 22 format (/'Energy-term weights (scaled):'//
720 & 'WSCC= ',f10.6,' (SC-SC)'/
721 & 'WSCP= ',f10.6,' (SC-p)'/
722 & 'WELEC= ',f10.6,' (p-p electr)'/
723 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
724 & 'WBOND= ',f10.6,' (stretching)'/
725 & 'WANG= ',f10.6,' (bending)'/
726 & 'WSCLOC= ',f10.6,' (SC local)'/
727 & 'WTOR= ',f10.6,' (torsional)'/
728 & 'WTORD= ',f10.6,' (double torsional)'/
729 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
730 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
731 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
732 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
733 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
734 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
735 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
736 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
737 & 'WTURN6= ',f10.6,' (turns, 6th order)')
738 if(me.eq.king.or..not.out1file)
739 & write (iout,*) "Reference temperature for weights calculation:",
741 call reada(weightcard,"D0CM",d0cm,3.78d0)
742 call reada(weightcard,"AKCM",akcm,15.1d0)
743 call reada(weightcard,"AKTH",akth,11.0d0)
744 call reada(weightcard,"AKCT",akct,12.0d0)
745 call reada(weightcard,"V1SS",v1ss,-1.08d0)
746 call reada(weightcard,"V2SS",v2ss,7.61d0)
747 call reada(weightcard,"V3SS",v3ss,13.7d0)
748 call reada(weightcard,"EBR",ebr,-5.50D0)
749 call reada(weightcard,"ATRISS",atriss,0.301D0)
750 call reada(weightcard,"BTRISS",btriss,0.021D0)
751 call reada(weightcard,"CTRISS",ctriss,1.001D0)
752 call reada(weightcard,"DTRISS",dtriss,1.001D0)
753 write (iout,*) "ATRISS=", atriss
754 write (iout,*) "BTRISS=", btriss
755 write (iout,*) "CTRISS=", ctriss
756 write (iout,*) "DTRISS=", dtriss
757 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
759 dyn_ss_mask(i)=.false.
763 dyn_ssbond_ij(i,j)=1.0d300
766 call reada(weightcard,"HT",Ht,0.0D0)
768 ss_depth=ebr/wsc-0.25*eps(1,1)
769 Ht=Ht/wsc-0.25*eps(1,1)
770 akcm=akcm*wstrain/wsc
771 akth=akth*wstrain/wsc
772 akct=akct*wstrain/wsc
773 v1ss=v1ss*wstrain/wsc
774 v2ss=v2ss*wstrain/wsc
775 v3ss=v3ss*wstrain/wsc
777 if (wstrain.ne.0.0) then
778 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
784 if(me.eq.king.or..not.out1file) then
785 write (iout,*) "Parameters of the SS-bond potential:"
786 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
788 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
789 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
790 write (iout,*)" HT",Ht
791 print *,'indpdb=',indpdb,' pdbref=',pdbref
793 if (indpdb.gt.0 .or. pdbref) then
794 read(inp,'(a)') pdbfile
795 if(me.eq.king.or..not.out1file)
796 & write (iout,'(2a)') 'PDB data will be read from file ',
797 & pdbfile(:ilen(pdbfile))
798 open(ipdbin,file=pdbfile,status='old',err=33)
800 33 write (iout,'(a)') 'Error opening PDB file.'
803 c write (iout,*) 'Begin reading pdb data'
806 c write (iout,*) 'Finished reading pdb data'
808 if(me.eq.king.or..not.out1file)
809 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
810 & ' nstart_sup=',nstart_sup
812 itype_pdb(i)=itype(i)
816 nct=nstart_sup+nsup-1
817 call contact(.false.,ncont_ref,icont_ref,co)
820 if(me.eq.king.or..not.out1file)
821 & write(iout,*)'Adding sidechains'
825 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
828 do while (fail.and.nsi.le.maxsi)
829 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
832 if(fail) write(iout,*)'Adding sidechain failed for res ',
833 & i,' after ',nsi,' trials'
838 if (indpdb.eq.0) then
839 C Read sequence if not taken from the pdb file.
841 c print *,'nres=',nres
842 if (iscode.gt.0) then
843 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
845 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
847 C Convert sequence to numeric code
849 itype(i)=rescode(i,sequence(i),iscode)
851 C Assign initial virtual bond lengths
857 vbld(i+nres)=dsc(iabs(itype(i)))
858 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
859 c write (iout,*) "i",i," itype",itype(i),
860 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
864 c print '(20i4)',(itype(i),i=1,nres)
867 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
869 if (itype(i).eq.ntyp1) then
873 else if (iabs(itype(i+1)).ne.20) then
875 else if (iabs(itype(i)).ne.20) then
882 if(me.eq.king.or..not.out1file)then
883 write (iout,*) "ITEL"
885 write (iout,*) i,itype(i),itel(i)
887 print *,'Call Read_Bridge.'
890 C 8/13/98 Set limits to generating the dihedral angles
895 read (inp,*) ndih_constr
896 if (ndih_constr.gt.0) then
898 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
900 if(me.eq.king.or..not.out1file)then
902 & 'There are',ndih_constr,' constraints on phi angles.'
904 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
909 phi0(i)=deg2rad*phi0(i)
910 drange(i)=deg2rad*drange(i)
912 C if(me.eq.king.or..not.out1file)
913 C & write (iout,*) 'FTORS',ftors
916 phibound(1,ii) = phi0(i)-drange(i)
917 phibound(2,ii) = phi0(i)+drange(i)
920 C first setting the theta boundaries to 0 to pi
921 C this mean that there is no energy penalty for any angle occuring this can be applied
922 C for generate random conformation but is not implemented in this way
927 C begin reading theta constrains this is quartic constrains allowing to
928 C have smooth second derivative
929 if (with_theta_constr) then
930 C with_theta_constr is keyword allowing for occurance of theta constrains
931 read (inp,*) ntheta_constr
932 C ntheta_constr is the number of theta constrains
933 if (ntheta_constr.gt.0) then
935 read (inp,*) (itheta_constr(i),theta_constr0(i),
936 & theta_drange(i),for_thet_constr(i),
938 C the above code reads from 1 to ntheta_constr
939 C itheta_constr(i) residue i for which is theta_constr
940 C theta_constr0 the global minimum value
941 C theta_drange is range for which there is no energy penalty
942 C for_thet_constr is the force constant for quartic energy penalty
944 if(me.eq.king.or..not.out1file)then
946 & 'There are',ntheta_constr,' constraints on phi angles.'
948 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
954 theta_constr0(i)=deg2rad*theta_constr0(i)
955 theta_drange(i)=deg2rad*theta_drange(i)
957 C if(me.eq.king.or..not.out1file)
958 C & write (iout,*) 'FTORS',ftors
959 C do i=1,ntheta_constr
960 C ii = itheta_constr(i)
961 C thetabound(1,ii) = phi0(i)-drange(i)
962 C thetabound(2,ii) = phi0(i)+drange(i)
964 endif ! ntheta_constr.gt.0
965 endif! with_theta_constr
967 C with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
968 C write (iout,*) "with_dihed_constr ",with_dihed_constr
973 write (iout,'(a)') 'Boundaries in phi angle sampling:'
975 write (iout,'(a3,i5,2f10.1)')
976 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
982 cd print *,'NNT=',NNT,' NCT=',NCT
983 if (itype(1).eq.ntyp1) nnt=2
984 if (itype(nres).eq.ntyp1) nct=nct-1
986 if(me.eq.king.or..not.out1file)
987 & write (iout,'(a,i3)') 'nsup=',nsup
989 if (nsup.le.(nct-nnt+1)) then
990 do i=0,nct-nnt+1-nsup
991 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
997 & 'Error - sequences to be superposed do not match.'
1000 do i=0,nsup-(nct-nnt+1)
1001 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1003 nstart_sup=nstart_sup+i
1009 & 'Error - sequences to be superposed do not match.'
1012 if (nsup.eq.0) nsup=nct-nnt
1013 if (nstart_sup.eq.0) nstart_sup=nnt
1014 if (nstart_seq.eq.0) nstart_seq=nnt
1015 if(me.eq.king.or..not.out1file)
1016 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1017 & ' nstart_seq=',nstart_seq
1019 c--- Zscore rms -------
1020 if (nz_start.eq.0) nz_start=nnt
1021 if (nz_end.eq.0 .and. nsup.gt.0) then
1023 else if (nz_end.eq.0) then
1026 if(me.eq.king.or..not.out1file)then
1027 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1028 write (iout,*) 'IZ_SC=',iz_sc
1030 c----------------------
1033 if (.not.pdbref) then
1034 call read_angles(inp,*38)
1036 38 write (iout,'(a)') 'Error reading reference structure.'
1038 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1039 stop 'Error reading reference structure'
1043 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1050 cref(j,i,kkk)=c(j,i)
1053 call contact(.true.,ncont_ref,icont_ref,co)
1057 write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1059 if (constr_dist.gt.0) call read_dist_constr
1060 write (iout,*) "After read_dist_constr nhpb",nhpb
1061 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1063 if(me.eq.king.or..not.out1file)
1064 & write (iout,*) 'Contact order:',co
1066 if(me.eq.king.or..not.out1file)
1067 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1070 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1072 if(me.eq.king.or..not.out1file)
1073 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1074 & icont_ref(1,i),' ',
1075 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1079 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1080 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1081 & modecalc.ne.10) then
1082 C If input structure hasn't been supplied from the PDB file read or generate
1084 if (iranconf.eq.0 .and. .not. extconf) then
1085 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1086 & write (iout,'(a)') 'Initial geometry will be read in.'
1088 read(inp,'(8f10.5)',end=36,err=36)
1089 & ((c(l,k),l=1,3),k=1,nres),
1090 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1091 write (iout,*) "Exit READ_CART"
1092 write (iout,'(8f10.5)')
1093 & ((c(l,k),l=1,3),k=1,nres),
1094 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1095 call int_from_cart1(.true.)
1096 write (iout,*) "Finish INT_TO_CART"
1099 dc(j,i)=c(j,i+1)-c(j,i)
1100 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1104 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1106 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1107 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1113 call read_angles(inp,*36)
1116 36 write (iout,'(a)') 'Error reading angle file.'
1118 call mpi_finalize( MPI_COMM_WORLD,IERR )
1120 stop 'Error reading angle file.'
1122 else if (extconf) then
1123 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1124 & write (iout,'(a)') 'Extended chain initial geometry.'
1126 theta(i)=90d0*deg2rad
1129 phi(i)=180d0*deg2rad
1132 alph(i)=110d0*deg2rad
1135 omeg(i)=-120d0*deg2rad
1136 if (itype(i).le.0) omeg(i)=-omeg(i)
1138 call chainbuild_extconf
1140 if(me.eq.king.or..not.out1file)
1141 & write (iout,'(a)') 'Random-generated initial geometry.'
1145 if (me.eq.king .or. fg_rank.eq.0 .and. (
1146 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1150 call gen_rand_conf(itmp,*30)
1152 30 write (iout,*) 'Failed to generate random conformation',
1153 & ', itrial=',itrial
1154 write (*,*) 'Processor:',me,
1155 & ' Failed to generate random conformation',
1165 write (iout,'(a,i3,a)') 'Processor:',me,
1166 & ' error in generating random conformation.'
1167 write (*,'(a,i3,a)') 'Processor:',me,
1168 & ' error in generating random conformation.'
1171 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1176 & ' error in generating random conformation.'
1181 elseif (modecalc.eq.4) then
1182 read (inp,'(a)') intinname
1183 open (intin,file=intinname,status='old',err=333)
1184 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1185 & write (iout,'(a)') 'intinname',intinname
1186 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1188 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1190 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1192 stop 'Error opening angle file.'
1196 C Generate distance constraints, if the PDB structure is to be regularized.
1197 if (nthread.gt.0) then
1198 call read_threadbase
1201 if (me.eq.king .or. .not. out1file)
1203 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1204 write (iout,'(/a,i3,a)')
1205 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1206 write (iout,'(20i4)') (iss(i),i=1,ns)
1208 write(iout,*)"Running with dynamic disulfide-bond formation"
1210 write (iout,'(/a/)') 'Pre-formed links are:'
1216 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1217 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1223 if (ns.gt.0.and.dyn_ss) then
1227 forcon(i-nss)=forcon(i)
1234 dyn_ss_mask(iss(i))=.true.
1237 if (i2ndstr.gt.0) call secstrp2dihc
1238 c call geom_to_var(nvar,x)
1239 c call etotal(energia(0))
1240 c call enerprint(energia(0))
1241 c call briefout(0,etot)
1243 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1244 cd write (iout,'(a)') 'Variable list:'
1245 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1247 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1248 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1249 & 'Processor',myrank,': end reading molecular data.'
1254 c--------------------------------------------------------------------------
1255 logical function seq_comp(itypea,itypeb,length)
1257 integer length,itypea(length),itypeb(length)
1260 if (itypea(i).ne.itypeb(i)) then
1268 c-----------------------------------------------------------------------------
1269 subroutine read_bridge
1270 C Read information about disulfide bridges.
1271 implicit real*8 (a-h,o-z)
1272 include 'DIMENSIONS'
1276 include 'COMMON.IOUNITS'
1277 include 'COMMON.GEO'
1278 include 'COMMON.VAR'
1279 include 'COMMON.INTERACT'
1280 include 'COMMON.LOCAL'
1281 include 'COMMON.NAMES'
1282 include 'COMMON.CHAIN'
1283 include 'COMMON.FFIELD'
1284 include 'COMMON.SBRIDGE'
1285 include 'COMMON.HEADER'
1286 include 'COMMON.CONTROL'
1287 include 'COMMON.DBASE'
1288 include 'COMMON.THREAD'
1289 include 'COMMON.TIME1'
1290 include 'COMMON.SETUP'
1291 C Read bridging residues.
1292 read (inp,*) ns,(iss(i),i=1,ns)
1294 if(me.eq.king.or..not.out1file)
1295 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1296 C Check whether the specified bridging residues are cystines.
1298 if (itype(iss(i)).ne.1) then
1299 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1300 & 'Do you REALLY think that the residue ',
1301 & restyp(itype(iss(i))),i,
1302 & ' can form a disulfide bridge?!!!'
1303 write (*,'(2a,i3,a)')
1304 & 'Do you REALLY think that the residue ',
1305 & restyp(itype(iss(i))),i,
1306 & ' can form a disulfide bridge?!!!'
1308 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1313 C Read preformed bridges.
1315 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1317 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1320 C Check if the residues involved in bridges are in the specified list of
1321 C bridging residues.
1324 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1325 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1326 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1327 & ' contains residues present in other pairs.'
1328 write (*,'(a,i3,a)') 'Disulfide pair',i,
1329 & ' contains residues present in other pairs.'
1331 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1337 if (ihpb(i).eq.iss(j)) goto 10
1339 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1342 if (jhpb(i).eq.iss(j)) goto 20
1344 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1350 ihpb(i)=ihpb(i)+nres
1351 jhpb(i)=jhpb(i)+nres
1357 c----------------------------------------------------------------------------
1358 subroutine read_x(kanal,*)
1359 implicit real*8 (a-h,o-z)
1360 include 'DIMENSIONS'
1361 include 'COMMON.GEO'
1362 include 'COMMON.VAR'
1363 include 'COMMON.CHAIN'
1364 include 'COMMON.IOUNITS'
1365 include 'COMMON.CONTROL'
1366 include 'COMMON.LOCAL'
1367 include 'COMMON.INTERACT'
1368 c Read coordinates from input
1370 read(kanal,'(8f10.5)',end=10,err=10)
1371 & ((c(l,k),l=1,3),k=1,nres),
1372 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1375 c(j,2*nres)=c(j,nres)
1377 call int_from_cart1(.false.)
1380 dc(j,i)=c(j,i+1)-c(j,i)
1381 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1385 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1387 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1388 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1396 c----------------------------------------------------------------------------
1397 subroutine read_threadbase
1398 implicit real*8 (a-h,o-z)
1399 include 'DIMENSIONS'
1400 include 'COMMON.IOUNITS'
1401 include 'COMMON.GEO'
1402 include 'COMMON.VAR'
1403 include 'COMMON.INTERACT'
1404 include 'COMMON.LOCAL'
1405 include 'COMMON.NAMES'
1406 include 'COMMON.CHAIN'
1407 include 'COMMON.FFIELD'
1408 include 'COMMON.SBRIDGE'
1409 include 'COMMON.HEADER'
1410 include 'COMMON.CONTROL'
1411 include 'COMMON.DBASE'
1412 include 'COMMON.THREAD'
1413 include 'COMMON.TIME1'
1414 C Read pattern database for threading.
1415 read (icbase,*) nseq
1417 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1418 & nres_base(2,i),nres_base(3,i)
1419 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1421 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1422 c & nres_base(2,i),nres_base(3,i)
1423 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1427 if (weidis.eq.0.0D0) weidis=0.1D0
1436 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1437 write (iout,'(a,i5)') 'nexcl: ',nexcl
1438 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1441 c------------------------------------------------------------------------------
1442 subroutine setup_var
1443 implicit real*8 (a-h,o-z)
1444 include 'DIMENSIONS'
1445 include 'COMMON.IOUNITS'
1446 include 'COMMON.GEO'
1447 include 'COMMON.VAR'
1448 include 'COMMON.INTERACT'
1449 include 'COMMON.LOCAL'
1450 include 'COMMON.NAMES'
1451 include 'COMMON.CHAIN'
1452 include 'COMMON.FFIELD'
1453 include 'COMMON.SBRIDGE'
1454 include 'COMMON.HEADER'
1455 include 'COMMON.CONTROL'
1456 include 'COMMON.DBASE'
1457 include 'COMMON.THREAD'
1458 include 'COMMON.TIME1'
1459 C Set up variable list.
1465 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1467 ialph(i,1)=nvar+nside
1471 if (indphi.gt.0) then
1473 else if (indback.gt.0) then
1478 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1481 c----------------------------------------------------------------------------
1482 subroutine gen_dist_constr
1483 C Generate CA distance constraints.
1484 implicit real*8 (a-h,o-z)
1485 include 'DIMENSIONS'
1486 include 'COMMON.IOUNITS'
1487 include 'COMMON.GEO'
1488 include 'COMMON.VAR'
1489 include 'COMMON.INTERACT'
1490 include 'COMMON.LOCAL'
1491 include 'COMMON.NAMES'
1492 include 'COMMON.CHAIN'
1493 include 'COMMON.FFIELD'
1494 include 'COMMON.SBRIDGE'
1495 include 'COMMON.HEADER'
1496 include 'COMMON.CONTROL'
1497 include 'COMMON.DBASE'
1498 include 'COMMON.THREAD'
1499 include 'COMMON.TIME1'
1500 dimension itype_pdb(maxres)
1501 common /pizda/ itype_pdb
1503 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1504 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1505 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1507 do i=nstart_sup,nstart_sup+nsup-1
1508 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1509 cd & ' seq_pdb', restyp(itype_pdb(i))
1510 do j=i+2,nstart_sup+nsup-1
1512 ihpb(nhpb)=i+nstart_seq-nstart_sup
1513 jhpb(nhpb)=j+nstart_seq-nstart_sup
1515 dhpb(nhpb)=dist(i,j)
1518 cd write (iout,'(a)') 'Distance constraints:'
1523 cd if (ii.gt.nres) then
1528 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1529 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1530 cd & dhpb(i),forcon(i)
1534 c----------------------------------------------------------------------------
1536 implicit real*8 (a-h,o-z)
1537 include 'DIMENSIONS'
1538 include 'COMMON.MAP'
1539 include 'COMMON.IOUNITS'
1540 character*3 angid(4) /'THE','PHI','ALP','OME'/
1541 character*80 mapcard,ucase
1543 read (inp,'(a)') mapcard
1544 mapcard=ucase(mapcard)
1545 if (index(mapcard,'PHI').gt.0) then
1547 else if (index(mapcard,'THE').gt.0) then
1549 else if (index(mapcard,'ALP').gt.0) then
1551 else if (index(mapcard,'OME').gt.0) then
1554 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1555 stop 'Error - illegal variable spec in MAP card.'
1557 call readi (mapcard,'RES1',res1(imap),0)
1558 call readi (mapcard,'RES2',res2(imap),0)
1559 if (res1(imap).eq.0) then
1560 res1(imap)=res2(imap)
1561 else if (res2(imap).eq.0) then
1562 res2(imap)=res1(imap)
1564 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1566 & 'Error - illegal definition of variable group in MAP.'
1567 stop 'Error - illegal definition of variable group in MAP.'
1569 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1570 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1571 call readi(mapcard,'NSTEP',nstep(imap),0)
1572 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1574 & 'Illegal boundary and/or step size specification in MAP.'
1575 stop 'Illegal boundary and/or step size specification in MAP.'
1580 c----------------------------------------------------------------------------
1582 implicit real*8 (a-h,o-z)
1583 include 'DIMENSIONS'
1584 include 'COMMON.IOUNITS'
1585 include 'COMMON.GEO'
1586 include 'COMMON.CSA'
1587 include 'COMMON.BANK'
1588 include 'COMMON.CONTROL'
1590 character*620 mcmcard
1591 call card_concat(mcmcard)
1593 call readi(mcmcard,'NCONF',nconf,50)
1594 call readi(mcmcard,'NADD',nadd,0)
1595 call readi(mcmcard,'JSTART',jstart,1)
1596 call readi(mcmcard,'JEND',jend,1)
1597 call readi(mcmcard,'NSTMAX',nstmax,500000)
1598 call readi(mcmcard,'N0',n0,1)
1599 call readi(mcmcard,'N1',n1,6)
1600 call readi(mcmcard,'N2',n2,4)
1601 call readi(mcmcard,'N3',n3,0)
1602 call readi(mcmcard,'N4',n4,0)
1603 call readi(mcmcard,'N5',n5,0)
1604 call readi(mcmcard,'N6',n6,10)
1605 call readi(mcmcard,'N7',n7,0)
1606 call readi(mcmcard,'N8',n8,0)
1607 call readi(mcmcard,'N9',n9,0)
1608 call readi(mcmcard,'N14',n14,0)
1609 call readi(mcmcard,'N15',n15,0)
1610 call readi(mcmcard,'N16',n16,0)
1611 call readi(mcmcard,'N17',n17,0)
1612 call readi(mcmcard,'N18',n18,0)
1614 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1616 call readi(mcmcard,'NDIFF',ndiff,2)
1617 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1618 call readi(mcmcard,'IS1',is1,1)
1619 call readi(mcmcard,'IS2',is2,8)
1620 call readi(mcmcard,'NRAN0',nran0,4)
1621 call readi(mcmcard,'NRAN1',nran1,2)
1622 call readi(mcmcard,'IRR',irr,1)
1623 call readi(mcmcard,'NSEED',nseed,20)
1624 call readi(mcmcard,'NTOTAL',ntotal,10000)
1625 call reada(mcmcard,'CUT1',cut1,2.0d0)
1626 call reada(mcmcard,'CUT2',cut2,5.0d0)
1627 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1628 call readi(mcmcard,'ICMAX',icmax,3)
1629 call readi(mcmcard,'IRESTART',irestart,0)
1630 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1633 call reada(mcmcard,'DELE',dele,20.0d0)
1634 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1635 call readi(mcmcard,'IREF',iref,0)
1636 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1637 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1638 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1639 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1640 write (iout,*) "NCONF_IN",nconf_in
1643 c----------------------------------------------------------------------------
1644 cfmc subroutine mcmfread
1645 cfmc implicit real*8 (a-h,o-z)
1646 cfmc include 'DIMENSIONS'
1647 cfmc include 'COMMON.MCMF'
1648 cfmc include 'COMMON.IOUNITS'
1649 cfmc include 'COMMON.GEO'
1650 cfmc character*80 ucase
1651 cfmc character*620 mcmcard
1652 cfmc call card_concat(mcmcard)
1654 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1655 cfmc write(iout,*)'MAXRANT=',maxrant
1656 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1657 cfmc write(iout,*)'MAXFAM=',maxfam
1658 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1659 cfmc write(iout,*)'NNET1=',nnet1
1660 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1661 cfmc write(iout,*)'NNET2=',nnet2
1662 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1663 cfmc write(iout,*)'NNET3=',nnet3
1664 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1665 cfmc write(iout,*)'ILASTT=',ilastt
1666 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1667 cfmc write(iout,*)'MAXSTR=',maxstr
1668 cfmc maxstr_f=maxstr/maxfam
1669 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1670 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1671 cfmc write(iout,*)'NMCMF=',nmcmf
1672 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1673 cfmc write(iout,*)'IFOCUS=',ifocus
1674 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1675 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1676 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1677 cfmc write(iout,*)'INTPRT=',intprt
1678 cfmc call readi(mcmcard,'IPRT',iprt,100)
1679 cfmc write(iout,*)'IPRT=',iprt
1680 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1681 cfmc write(iout,*)'IMAXTR=',imaxtr
1682 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1683 cfmc write(iout,*)'MAXEVEN=',maxeven
1684 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1685 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1686 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1687 cfmc write(iout,*)'INIMIN=',inimin
1688 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1689 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1690 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1691 cfmc write(iout,*)'NTHREAD=',nthread
1692 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1693 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1694 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1695 cfmc write(iout,*)'MAXPERT=',maxpert
1696 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1697 cfmc write(iout,*)'IRMSD=',irmsd
1698 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1699 cfmc write(iout,*)'DENEMIN=',denemin
1700 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1701 cfmc write(iout,*)'RCUT1S=',rcut1s
1702 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1703 cfmc write(iout,*)'RCUT1E=',rcut1e
1704 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1705 cfmc write(iout,*)'RCUT2S=',rcut2s
1706 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1707 cfmc write(iout,*)'RCUT2E=',rcut2e
1708 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1709 cfmc write(iout,*)'DPERT1=',d_pert1
1710 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1711 cfmc write(iout,*)'DPERT1A=',d_pert1a
1712 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1713 cfmc write(iout,*)'DPERT2=',d_pert2
1714 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1715 cfmc write(iout,*)'DPERT2A=',d_pert2a
1716 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1717 cfmc write(iout,*)'DPERT2B=',d_pert2b
1718 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1719 cfmc write(iout,*)'DPERT2C=',d_pert2c
1720 cfmc d_pert1=deg2rad*d_pert1
1721 cfmc d_pert1a=deg2rad*d_pert1a
1722 cfmc d_pert2=deg2rad*d_pert2
1723 cfmc d_pert2a=deg2rad*d_pert2a
1724 cfmc d_pert2b=deg2rad*d_pert2b
1725 cfmc d_pert2c=deg2rad*d_pert2c
1726 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1727 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1728 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1729 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1730 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1731 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1732 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1733 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1734 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1735 cfmc write(iout,*)'RCUTINI=',rcutini
1736 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1737 cfmc write(iout,*)'GRAT=',grat
1738 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1739 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1743 c----------------------------------------------------------------------------
1745 implicit real*8 (a-h,o-z)
1746 include 'DIMENSIONS'
1747 include 'COMMON.MCM'
1748 include 'COMMON.MCE'
1749 include 'COMMON.IOUNITS'
1751 character*320 mcmcard
1752 call card_concat(mcmcard)
1753 call readi(mcmcard,'MAXACC',maxacc,100)
1754 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1755 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1756 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1757 call readi(mcmcard,'MAXREPM',maxrepm,200)
1758 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1759 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1760 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1761 call reada(mcmcard,'E_UP',e_up,5.0D0)
1762 call reada(mcmcard,'DELTE',delte,0.1D0)
1763 call readi(mcmcard,'NSWEEP',nsweep,5)
1764 call readi(mcmcard,'NSTEPH',nsteph,0)
1765 call readi(mcmcard,'NSTEPC',nstepc,0)
1766 call reada(mcmcard,'TMIN',tmin,298.0D0)
1767 call reada(mcmcard,'TMAX',tmax,298.0D0)
1768 call readi(mcmcard,'NWINDOW',nwindow,0)
1769 call readi(mcmcard,'PRINT_MC',print_mc,0)
1770 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1771 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1772 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1773 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1774 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1775 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1776 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1777 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1778 if (nwindow.gt.0) then
1779 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1781 winlen(i)=winend(i)-winstart(i)+1
1784 if (tmax.lt.tmin) tmax=tmin
1785 if (tmax.eq.tmin) then
1789 if (nstepc.gt.0 .and. nsteph.gt.0) then
1790 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1791 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1793 C Probabilities of different move types
1794 sumpro_type(0)=0.0D0
1795 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1796 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1797 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1798 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1799 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1800 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1801 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1803 print *,'i',i,' sumprotype',sumpro_type(i)
1804 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1805 print *,'i',i,' sumprotype',sumpro_type(i)
1809 c----------------------------------------------------------------------------
1810 subroutine read_minim
1811 implicit real*8 (a-h,o-z)
1812 include 'DIMENSIONS'
1813 include 'COMMON.MINIM'
1814 include 'COMMON.IOUNITS'
1816 character*320 minimcard
1817 call card_concat(minimcard)
1818 call readi(minimcard,'MAXMIN',maxmin,2000)
1819 call readi(minimcard,'MAXFUN',maxfun,5000)
1820 call readi(minimcard,'MINMIN',minmin,maxmin)
1821 call readi(minimcard,'MINFUN',minfun,maxmin)
1822 call reada(minimcard,'TOLF',tolf,1.0D-2)
1823 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1824 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1825 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1826 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1827 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1828 & 'Options in energy minimization:'
1829 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1830 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1831 & 'MinMin:',MinMin,' MinFun:',MinFun,
1832 & ' TolF:',TolF,' RTolF:',RTolF
1835 c----------------------------------------------------------------------------
1836 subroutine read_angles(kanal,*)
1837 implicit real*8 (a-h,o-z)
1838 include 'DIMENSIONS'
1839 include 'COMMON.GEO'
1840 include 'COMMON.VAR'
1841 include 'COMMON.CHAIN'
1842 include 'COMMON.IOUNITS'
1843 include 'COMMON.CONTROL'
1844 c Read angles from input
1846 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1847 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1848 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1849 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1852 c 9/7/01 avoid 180 deg valence angle
1853 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1855 theta(i)=deg2rad*theta(i)
1856 phi(i)=deg2rad*phi(i)
1857 alph(i)=deg2rad*alph(i)
1858 omeg(i)=deg2rad*omeg(i)
1863 c----------------------------------------------------------------------------
1864 subroutine reada(rekord,lancuch,wartosc,default)
1866 character*(*) rekord,lancuch
1867 double precision wartosc,default
1870 iread=index(rekord,lancuch)
1871 if (iread.eq.0) then
1875 iread=iread+ilen(lancuch)+1
1876 read (rekord(iread:),*,err=10,end=10) wartosc
1881 c----------------------------------------------------------------------------
1882 subroutine readi(rekord,lancuch,wartosc,default)
1884 character*(*) rekord,lancuch
1885 integer wartosc,default
1888 iread=index(rekord,lancuch)
1889 if (iread.eq.0) then
1893 iread=iread+ilen(lancuch)+1
1894 read (rekord(iread:),*,err=10,end=10) wartosc
1899 c----------------------------------------------------------------------------
1900 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1903 integer tablica(dim),default
1904 character*(*) rekord,lancuch
1911 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1912 if (iread.eq.0) return
1913 iread=iread+ilen(lancuch)+1
1914 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1917 c----------------------------------------------------------------------------
1918 subroutine multreada(rekord,lancuch,tablica,dim,default)
1921 double precision tablica(dim),default
1922 character*(*) rekord,lancuch
1929 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1930 if (iread.eq.0) return
1931 iread=iread+ilen(lancuch)+1
1932 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1935 c----------------------------------------------------------------------------
1936 subroutine openunits
1937 implicit real*8 (a-h,o-z)
1938 include 'DIMENSIONS'
1941 character*16 form,nodename
1944 include 'COMMON.SETUP'
1945 include 'COMMON.IOUNITS'
1947 include 'COMMON.CONTROL'
1948 integer lenpre,lenpot,ilen,lentmp
1950 character*3 out1file_text,ucase
1953 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1954 call getenv_loc("PREFIX",prefix)
1956 call getenv_loc("POT",pot)
1957 call getenv_loc("DIRTMP",tmpdir)
1958 call getenv_loc("CURDIR",curdir)
1959 call getenv_loc("OUT1FILE",out1file_text)
1960 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1961 out1file_text=ucase(out1file_text)
1962 if (out1file_text(1:1).eq."Y") then
1965 out1file=fg_rank.gt.0
1970 if (lentmp.gt.0) then
1971 write (*,'(80(1h!))')
1972 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1973 write (*,'(80(1h!))')
1974 write (*,*)"All output files will be on node /tmp directory."
1976 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1977 if (me.eq.king) then
1978 write (*,*) "The master node is ",nodename
1979 else if (fg_rank.eq.0) then
1980 write (*,*) "I am the CG slave node ",nodename
1982 write (*,*) "I am the FG slave node ",nodename
1985 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1986 lenpre = lentmp+lenpre+1
1988 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1989 C Get the names and open the input files
1990 #if defined(WINIFL) || defined(WINPGI)
1991 open(1,file=pref_orig(:ilen(pref_orig))//
1992 & '.inp',status='old',readonly,shared)
1993 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1994 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1995 C Get parameter filenames and open the parameter files.
1996 call getenv_loc('BONDPAR',bondname)
1997 open (ibond,file=bondname,status='old',readonly,shared)
1998 call getenv_loc('THETPAR',thetname)
1999 open (ithep,file=thetname,status='old',readonly,shared)
2000 call getenv_loc('ROTPAR',rotname)
2001 open (irotam,file=rotname,status='old',readonly,shared)
2002 call getenv_loc('TORPAR',torname)
2003 open (itorp,file=torname,status='old',readonly,shared)
2004 call getenv_loc('TORDPAR',tordname)
2005 open (itordp,file=tordname,status='old',readonly,shared)
2006 call getenv_loc('TORKCC',torkccname)
2007 open (itorkcc,file=torkccname,status='old',readonly,shared)
2008 call getenv_loc('THETKCC',thetkccname)
2009 open (ithetkcc,file=thetkccname,status='old',readonly,shared)
2010 call getenv_loc('FOURIER',fouriername)
2011 open (ifourier,file=fouriername,status='old',readonly,shared)
2012 call getenv_loc('ELEPAR',elename)
2013 open (ielep,file=elename,status='old',readonly,shared)
2014 call getenv_loc('SIDEPAR',sidename)
2015 open (isidep,file=sidename,status='old',readonly,shared)
2016 call getenv_loc('LIPTRANPAR',liptranname)
2017 open (iliptranpar,file=liptranname,status='old',readonly,shared)
2018 #elif (defined CRAY) || (defined AIX)
2019 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2021 c print *,"Processor",myrank," opened file 1"
2022 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2023 c print *,"Processor",myrank," opened file 9"
2024 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2025 C Get parameter filenames and open the parameter files.
2026 call getenv_loc('BONDPAR',bondname)
2027 open (ibond,file=bondname,status='old',action='read')
2028 c print *,"Processor",myrank," opened file IBOND"
2029 call getenv_loc('THETPAR',thetname)
2030 open (ithep,file=thetname,status='old',action='read')
2031 c print *,"Processor",myrank," opened file ITHEP"
2032 call getenv_loc('ROTPAR',rotname)
2033 open (irotam,file=rotname,status='old',action='read')
2034 c print *,"Processor",myrank," opened file IROTAM"
2035 call getenv_loc('TORPAR',torname)
2036 open (itorp,file=torname,status='old',action='read')
2037 c print *,"Processor",myrank," opened file ITORP"
2038 call getenv_loc('TORDPAR',tordname)
2039 open (itordp,file=tordname,status='old',action='read')
2040 call getenv_loc('TORKCC',torkccname)
2041 open (itorkcc,file=torkccname,status='old',action='read')
2042 call getenv_loc('THETKCC',thetkccname)
2043 open (ithetkcc,file=thetkccname,status='old',action='read')
2044 c print *,"Processor",myrank," opened file ITORDP"
2045 call getenv_loc('SCCORPAR',sccorname)
2046 open (isccor,file=sccorname,status='old',action='read')
2047 c print *,"Processor",myrank," opened file ISCCOR"
2048 call getenv_loc('FOURIER',fouriername)
2049 open (ifourier,file=fouriername,status='old',action='read')
2050 c print *,"Processor",myrank," opened file IFOURIER"
2051 call getenv_loc('ELEPAR',elename)
2052 open (ielep,file=elename,status='old',action='read')
2053 c print *,"Processor",myrank," opened file IELEP"
2054 call getenv_loc('SIDEPAR',sidename)
2055 open (isidep,file=sidename,status='old',action='read')
2056 call getenv_loc('LIPTRANPAR',liptranname)
2057 open (iliptranpar,file=liptranname,status='old',action='read')
2058 c print *,"Processor",myrank," opened file ISIDEP"
2059 c print *,"Processor",myrank," opened parameter files"
2061 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2062 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2063 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2064 C Get parameter filenames and open the parameter files.
2065 call getenv_loc('BONDPAR',bondname)
2066 open (ibond,file=bondname,status='old')
2067 call getenv_loc('THETPAR',thetname)
2068 open (ithep,file=thetname,status='old')
2069 call getenv_loc('ROTPAR',rotname)
2070 open (irotam,file=rotname,status='old')
2071 call getenv_loc('TORPAR',torname)
2072 open (itorp,file=torname,status='old')
2073 call getenv_loc('TORDPAR',tordname)
2074 open (itordp,file=tordname,status='old')
2075 call getenv_loc('TORKCC',torkccname)
2076 open (itorkcc,file=torkccname,status='old')
2077 call getenv_loc('THETKCC',thetkccname)
2078 open (ithetkcc,file=thetkccname,status='old')
2079 call getenv_loc('SCCORPAR',sccorname)
2080 open (isccor,file=sccorname,status='old')
2081 call getenv_loc('FOURIER',fouriername)
2082 open (ifourier,file=fouriername,status='old')
2083 call getenv_loc('ELEPAR',elename)
2084 open (ielep,file=elename,status='old')
2085 call getenv_loc('SIDEPAR',sidename)
2086 open (isidep,file=sidename,status='old')
2087 call getenv_loc('LIPTRANPAR',liptranname)
2088 open (iliptranpar,file=liptranname,status='old')
2090 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2092 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2093 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2094 C Get parameter filenames and open the parameter files.
2095 call getenv_loc('BONDPAR',bondname)
2096 open (ibond,file=bondname,status='old',readonly)
2097 call getenv_loc('THETPAR',thetname)
2098 open (ithep,file=thetname,status='old',readonly)
2099 call getenv_loc('ROTPAR',rotname)
2100 open (irotam,file=rotname,status='old',readonly)
2101 call getenv_loc('TORPAR',torname)
2102 open (itorp,file=torname,status='old',readonly)
2103 call getenv_loc('TORDPAR',tordname)
2104 open (itordp,file=tordname,status='old',readonly)
2105 call getenv_loc('TORKCC',torkccname)
2106 open (itorkcc,file=torkccname,status='old',readonly)
2107 call getenv_loc('THETKCC',thetkccname)
2108 open (ithetkcc,file=thetkccname,status='old',readonly)
2109 call getenv_loc('SCCORPAR',sccorname)
2110 open (isccor,file=sccorname,status='old',readonly)
2112 call getenv_loc('THETPARPDB',thetname_pdb)
2113 print *,"thetname_pdb ",thetname_pdb
2114 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2115 print *,ithep_pdb," opened"
2117 call getenv_loc('FOURIER',fouriername)
2118 open (ifourier,file=fouriername,status='old',readonly)
2119 call getenv_loc('ELEPAR',elename)
2120 open (ielep,file=elename,status='old',readonly)
2121 call getenv_loc('SIDEPAR',sidename)
2122 open (isidep,file=sidename,status='old',readonly)
2123 call getenv_loc('LIPTRANPAR',liptranname)
2124 open (iliptranpar,file=liptranname,status='old',action='read')
2126 call getenv_loc('ROTPARPDB',rotname_pdb)
2127 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2132 C 8/9/01 In the newest version SCp interaction constants are read from a file
2133 C Use -DOLDSCP to use hard-coded constants instead.
2135 call getenv_loc('SCPPAR',scpname)
2136 #if defined(WINIFL) || defined(WINPGI)
2137 open (iscpp,file=scpname,status='old',readonly,shared)
2138 #elif (defined CRAY) || (defined AIX)
2139 open (iscpp,file=scpname,status='old',action='read')
2141 open (iscpp,file=scpname,status='old')
2143 open (iscpp,file=scpname,status='old',readonly)
2146 call getenv_loc('PATTERN',patname)
2147 #if defined(WINIFL) || defined(WINPGI)
2148 open (icbase,file=patname,status='old',readonly,shared)
2149 #elif (defined CRAY) || (defined AIX)
2150 open (icbase,file=patname,status='old',action='read')
2152 open (icbase,file=patname,status='old')
2154 open (icbase,file=patname,status='old',readonly)
2157 C Open output file only for CG processes
2158 c print *,"Processor",myrank," fg_rank",fg_rank
2159 if (fg_rank.eq.0) then
2161 if (nodes.eq.1) then
2164 npos = dlog10(dfloat(nodes-1))+1
2166 if (npos.lt.3) npos=3
2167 write (liczba,'(i1)') npos
2168 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2170 write (liczba,form) me
2171 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2172 & liczba(:ilen(liczba))
2173 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2175 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2177 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2178 & liczba(:ilen(liczba))//'.mol2'
2179 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2180 & liczba(:ilen(liczba))//'.stat'
2182 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2183 & //liczba(:ilen(liczba))//'.stat')
2184 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2187 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2188 & liczba(:ilen(liczba))//'.const'
2193 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2194 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2195 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2196 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2197 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2199 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2201 rest2name=prefix(:ilen(prefix))//'.rst'
2203 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2206 #if defined(AIX) || defined(PGI)
2207 if (me.eq.king .or. .not. out1file)
2208 & open(iout,file=outname,status='unknown')
2210 if (fg_rank.gt.0) then
2211 write (liczba,'(i3.3)') myrank/nfgtasks
2212 write (ll,'(bz,i3.3)') fg_rank
2213 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2218 open(igeom,file=intname,status='unknown',position='append')
2219 open(ipdb,file=pdbname,status='unknown')
2220 open(imol2,file=mol2name,status='unknown')
2221 open(istat,file=statname,status='unknown',position='append')
2223 c1out open(iout,file=outname,status='unknown')
2226 if (me.eq.king .or. .not.out1file)
2227 & open(iout,file=outname,status='unknown')
2229 if (fg_rank.gt.0) then
2230 write (liczba,'(i3.3)') myrank/nfgtasks
2231 write (ll,'(bz,i3.3)') fg_rank
2232 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2237 open(igeom,file=intname,status='unknown',access='append')
2238 open(ipdb,file=pdbname,status='unknown')
2239 open(imol2,file=mol2name,status='unknown')
2240 open(istat,file=statname,status='unknown',access='append')
2242 c1out open(iout,file=outname,status='unknown')
2245 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2246 csa_seed=prefix(:lenpre)//'.CSA.seed'
2247 csa_history=prefix(:lenpre)//'.CSA.history'
2248 csa_bank=prefix(:lenpre)//'.CSA.bank'
2249 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2250 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2251 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2252 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2253 csa_int=prefix(:lenpre)//'.int'
2254 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2255 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2256 csa_in=prefix(:lenpre)//'.CSA.in'
2257 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2260 write (iout,'(80(1h-))')
2261 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2262 write (iout,'(80(1h-))')
2263 write (iout,*) "Input file : ",
2264 & pref_orig(:ilen(pref_orig))//'.inp'
2265 write (iout,*) "Output file : ",
2266 & outname(:ilen(outname))
2268 write (iout,*) "Sidechain potential file : ",
2269 & sidename(:ilen(sidename))
2271 write (iout,*) "SCp potential file : ",
2272 & scpname(:ilen(scpname))
2274 write (iout,*) "Electrostatic potential file : ",
2275 & elename(:ilen(elename))
2276 write (iout,*) "Cumulant coefficient file : ",
2277 & fouriername(:ilen(fouriername))
2278 write (iout,*) "Torsional parameter file : ",
2279 & torname(:ilen(torname))
2280 write (iout,*) "Double torsional parameter file : ",
2281 & tordname(:ilen(tordname))
2282 write (iout,*) "SCCOR parameter file : ",
2283 & sccorname(:ilen(sccorname))
2284 write (iout,*) "Bond & inertia constant file : ",
2285 & bondname(:ilen(bondname))
2286 write (iout,*) "Bending parameter file : ",
2287 & thetname(:ilen(thetname))
2288 write (iout,*) "Rotamer parameter file : ",
2289 & rotname(:ilen(rotname))
2290 write (iout,*) "Thetpdb parameter file : ",
2291 & thetname_pdb(:ilen(thetname_pdb))
2292 write (iout,*) "Threading database : ",
2293 & patname(:ilen(patname))
2295 &write (iout,*)" DIRTMP : ",
2297 write (iout,'(80(1h-))')
2301 c----------------------------------------------------------------------------
2302 subroutine card_concat(card)
2303 implicit real*8 (a-h,o-z)
2304 include 'DIMENSIONS'
2305 include 'COMMON.IOUNITS'
2307 character*80 karta,ucase
2309 read (inp,'(a)') karta
2312 do while (karta(80:80).eq.'&')
2313 card=card(:ilen(card)+1)//karta(:79)
2314 read (inp,'(a)') karta
2317 card=card(:ilen(card)+1)//karta
2320 c----------------------------------------------------------------------------------
2322 implicit real*8 (a-h,o-z)
2323 include 'DIMENSIONS'
2324 include 'COMMON.CHAIN'
2325 include 'COMMON.IOUNITS'
2327 open(irest2,file=rest2name,status='unknown')
2328 read(irest2,*) totT,EK,potE,totE,t_bath
2331 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2334 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2337 read (irest2,*) iset
2342 c---------------------------------------------------------------------------------
2343 subroutine read_fragments
2344 implicit real*8 (a-h,o-z)
2345 include 'DIMENSIONS'
2349 include 'COMMON.SETUP'
2350 include 'COMMON.CHAIN'
2351 include 'COMMON.IOUNITS'
2353 include 'COMMON.CONTROL'
2354 read(inp,*) nset,nfrag,npair,nfrag_back
2355 if(me.eq.king.or..not.out1file)
2356 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2357 & " nfrag_back",nfrag_back
2359 read(inp,*) mset(iset)
2361 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2363 if(me.eq.king.or..not.out1file)
2364 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2365 & ifrag(2,i,iset), qinfrag(i,iset)
2368 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2370 if(me.eq.king.or..not.out1file)
2371 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2372 & ipair(2,i,iset), qinpair(i,iset)
2375 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2376 & wfrag_back(3,i,iset),
2377 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2378 if(me.eq.king.or..not.out1file)
2379 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2380 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2385 C---------------------------------------------------------------------------
2386 subroutine read_afminp
2387 implicit real*8 (a-h,o-z)
2388 include 'DIMENSIONS'
2392 include 'COMMON.SETUP'
2393 include 'COMMON.CONTROL'
2394 include 'COMMON.CHAIN'
2395 include 'COMMON.IOUNITS'
2396 include 'COMMON.SBRIDGE'
2397 character*320 afmcard
2399 call card_concat(afmcard)
2400 call readi(afmcard,"BEG",afmbeg,0)
2401 call readi(afmcard,"END",afmend,0)
2402 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2403 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2404 print *,'FORCE=' ,forceAFMconst
2405 CCCC NOW PROPERTIES FOR AFM
2408 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2410 distafminit=dsqrt(distafminit)
2411 print *,'initdist',distafminit
2415 c-------------------------------------------------------------------------------
2416 subroutine read_dist_constr
2417 implicit real*8 (a-h,o-z)
2418 include 'DIMENSIONS'
2422 include 'COMMON.SETUP'
2423 include 'COMMON.CONTROL'
2424 include 'COMMON.CHAIN'
2425 include 'COMMON.IOUNITS'
2426 include 'COMMON.SBRIDGE'
2427 integer ifrag_(2,100),ipair_(2,100)
2428 double precision wfrag_(100),wpair_(100)
2429 character*500 controlcard
2431 write (iout,*) "Calling read_dist_constr"
2432 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2434 call card_concat(controlcard)
2435 call readi(controlcard,"NFRAG",nfrag_,0)
2436 call readi(controlcard,"NPAIR",npair_,0)
2437 call readi(controlcard,"NDIST",ndist_,0)
2438 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2439 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2440 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2441 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2442 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2443 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2444 c write (iout,*) "IFRAG"
2446 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2448 c write (iout,*) "IPAIR"
2450 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2454 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2455 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2456 & ifrag_(2,i)=nstart_sup+nsup-1
2457 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2459 if (wfrag_(i).gt.0.0d0) then
2460 do j=ifrag_(1,i),ifrag_(2,i)-1
2461 do k=j+1,ifrag_(2,i)
2462 c write (iout,*) "j",j," k",k
2464 if (constr_dist.eq.1) then
2469 forcon(nhpb)=wfrag_(i)
2470 else if (constr_dist.eq.2) then
2471 if (ddjk.le.dist_cut) then
2476 forcon(nhpb)=wfrag_(i)
2483 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2486 if (.not.out1file .or. me.eq.king)
2487 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2488 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2490 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2491 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2498 if (wpair_(i).gt.0.0d0) then
2506 do j=ifrag_(1,ii),ifrag_(2,ii)
2507 do k=ifrag_(1,jj),ifrag_(2,jj)
2511 forcon(nhpb)=wpair_(i)
2512 dhpb(nhpb)=dist(j,k)
2514 if (.not.out1file .or. me.eq.king)
2515 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2516 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2518 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2519 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2527 if (constr_dist.eq.11) then
2528 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2529 & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
2530 fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2533 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2534 & ibecarb(i),forcon(nhpb+1)
2536 if (forcon(nhpb+1).gt.0.0d0) then
2538 if (ibecarb(i).gt.0) then
2539 ihpb(i)=ihpb(i)+nres
2540 jhpb(i)=jhpb(i)+nres
2542 if (dhpb(nhpb).eq.0.0d0)
2543 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2545 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2546 C if (forcon(nhpb+1).gt.0.0d0) then
2548 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2550 if (.not.out1file .or. me.eq.king)
2551 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2552 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2554 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2555 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2562 c-------------------------------------------------------------------------------
2564 subroutine flush(iu)
2569 subroutine flush(iu)
2574 c------------------------------------------------------------------------------
2575 subroutine copy_to_tmp(source)
2576 include "DIMENSIONS"
2577 include "COMMON.IOUNITS"
2578 character*(*) source
2579 character* 256 tmpfile
2583 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2584 inquire(file=tmpfile,exist=ex)
2586 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2587 & " to temporary directory..."
2588 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2589 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2593 c------------------------------------------------------------------------------
2594 subroutine move_from_tmp(source)
2595 include "DIMENSIONS"
2596 include "COMMON.IOUNITS"
2597 character*(*) source
2600 write (*,*) "Moving ",source(:ilen(source)),
2601 & " from temporary directory to working directory"
2602 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2603 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2606 c------------------------------------------------------------------------------
2607 subroutine random_init(seed)
2609 C Initialize random number generator
2611 implicit real*8 (a-h,o-z)
2612 include 'DIMENSIONS'
2615 logical OKRandom, prng_restart
2617 integer iseed_array(4)
2619 include 'COMMON.IOUNITS'
2620 include 'COMMON.TIME1'
2621 include 'COMMON.THREAD'
2622 include 'COMMON.SBRIDGE'
2623 include 'COMMON.CONTROL'
2624 include 'COMMON.MCM'
2625 include 'COMMON.MAP'
2626 include 'COMMON.HEADER'
2627 include 'COMMON.CSA'
2628 include 'COMMON.CHAIN'
2629 include 'COMMON.MUCA'
2631 include 'COMMON.FFIELD'
2632 include 'COMMON.SETUP'
2633 iseed=-dint(dabs(seed))
2634 if (iseed.eq.0) then
2635 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2636 & 'Random seed undefined. The program will stop.'
2637 write (*,'(/80(1h*)/20x,a/80(1h*))')
2638 & 'Random seed undefined. The program will stop.'
2640 call mpi_finalize(mpi_comm_world,ierr)
2642 stop 'Bad random seed.'
2645 if (fg_rank.eq.0) then
2649 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2650 OKRandom = prng_restart(me,iseed)
2653 tmp=65536.0d0**(4-i)
2654 iseed_array(i) = dint(seed/tmp)
2655 seed=seed-iseed_array(i)*tmp
2658 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2659 & (iseed_array(i),i=1,4)
2660 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2661 & (iseed_array(i),i=1,4)
2662 OKRandom = prng_restart(me,iseed_array)
2665 c r1 = prng_next(me)
2666 r1=ran_number(0.0D0,1.0D0)
2668 & write (iout,*) 'ran_num',r1
2669 if (r1.lt.0.0d0) OKRandom=.false.
2671 if (.not.OKRandom) then
2672 write (iout,*) 'PRNG IS NOT WORKING!!!'
2673 print *,'PRNG IS NOT WORKING!!!'
2676 call mpi_abort(mpi_comm_world,error_msg,ierr)
2679 write (iout,*) 'too many processors for parallel prng'
2680 write (*,*) 'too many processors for parallel prng'
2688 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)