8 use control_data, only:maxterm_sccor
10 !-----------------------------------------------------------------------------
11 ! Max. number of residue types and parameters in expressions for
12 ! virtual-bond angle bending potentials
13 ! integer,parameter :: maxthetyp=3
14 ! integer,parameter :: maxthetyp1=maxthetyp+1
16 ! & maxtheterm2=6,maxtheterm3=4,maxsingle=6,maxdouble=4,
17 ! & mmaxtheterm=maxtheterm)
18 !-----------------------------------------------------------------------------
19 ! Max. number of types of dihedral angles & multiplicity of torsional barriers
20 ! and the number of terms in double torsionals
21 ! integer,parameter :: maxlor=3,maxtermd_1=8,maxtermd_2=8
22 ! parameter (maxtor=4,maxterm=10)
23 !-----------------------------------------------------------------------------
24 ! Max number of torsional terms in SCCOR
25 !el integer,parameter :: maxterm_sccor=6
26 !-----------------------------------------------------------------------------
27 character(len=1),dimension(:),allocatable :: secstruc !(maxres)
28 !-----------------------------------------------------------------------------
31 !-----------------------------------------------------------------------------
33 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
34 !-----------------------------------------------------------------------------
36 !-----------------------------------------------------------------------------
37 subroutine write_rbank(jlee,adif,nft)
40 use geometry_data, only: nres,rad2deg
41 ! implicit real*8 (a-h,o-z)
42 ! include 'DIMENSIONS'
43 ! include 'COMMON.IOUNITS'
44 ! include 'COMMON.CSA'
45 ! include 'COMMON.BANK'
46 ! include 'COMMON.CHAIN'
47 ! include 'COMMON.GEO'
49 integer :: nft,i,k,j,l,jlee
52 open(icsa_rbank,file=csa_rbank,status="unknown")
53 write (icsa_rbank,900) jlee,nbank,nstep,nft,icycle,adif
55 write (icsa_rbank,952) k,rene(k),rrmsn(k),rpncn(k)
58 write (icsa_rbank,850) (rad2deg*rvar(i,l,j,k),i=1,4)
65 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
67 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
71 end subroutine write_rbank
72 !-----------------------------------------------------------------------------
73 subroutine read_rbank(jlee,adif)
76 use geometry_data, only: nres,deg2rad
78 ! implicit real*8 (a-h,o-z)
79 ! include 'DIMENSIONS'
81 ! include 'COMMON.IOUNITS'
82 ! include 'COMMON.CSA'
83 ! include 'COMMON.BANK'
84 ! include 'COMMON.CHAIN'
85 ! include 'COMMON.GEO'
86 ! include 'COMMON.SETUP'
87 character(len=80) :: karta
89 integer :: nbankr,nstepr,nftr,icycler,kk,k,j,l,i,&
90 ierror,ierrcode,jlee,jleer
93 open(icsa_rbank,file=csa_rbank,status="old")
94 read (icsa_rbank,901) jleer,nbankr,nstepr,nftr,icycler,adif
95 print *,jleer,nbankr,nstepr,nftr,icycler,adif
96 ! print *, 'adif from read_rbank ',adif
98 if(nbankr.ne.nbank) then
99 write (iout,*) 'ERROR in READ_BANK: NBANKR',nbankr,&
101 call mpi_abort(mpi_comm_world,ierror,ierrcode)
103 if(jleer.ne.jlee) then
104 write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,&
106 call mpi_abort(mpi_comm_world,ierror,ierrcode)
112 read (icsa_rbank,'(a80)') karta
113 write(iout,*) "READ_RBANK: kk=",kk
115 ! if (index(karta,"*").gt.0) then
116 ! write (iout,*) "***** Stars in bankr ***** k=",k,
120 ! read (30,850) (rdummy,i=1,4)
125 call reada(karta,"total E",rene(kk),1.0d20)
126 call reada(karta,"rmsd from N",rrmsn(kk),0.0d0)
127 call reada(karta,"%NC",rpncn(kk),0.0d0)
128 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
129 "%NC",bpncn(kk),ibank(kk)
130 ! read (icsa_rbank,953) kdummy,rene(kk),rrmsn(kk),rpncn(kk)
133 read (icsa_rbank,850) (rvar(i,l,j,kk),i=1,4)
134 ! write (iout,850) (rvar(i,l,j,kk),i=1,4)
136 rvar(i,l,j,kk)=deg2rad*rvar(i,l,j,kk)
142 !d write (*,*) "read_rbank ******************* kk",kk,
144 if (kk.lt.nbankr) nbankr=kk
149 !d write (*,850) (rvar(i,l,j,kk),i=1,4)
156 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
157 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2)
160 end subroutine read_rbank
161 !-----------------------------------------------------------------------------
162 subroutine write_bank(jlee,nft)
165 use control_data, only: vdisulf
166 use geometry_data, only: nres,rad2deg
167 ! implicit real*8 (a-h,o-z)
168 ! include 'DIMENSIONS'
169 ! include 'COMMON.IOUNITS'
170 ! include 'COMMON.CSA'
171 ! include 'COMMON.BANK'
172 ! include 'COMMON.CHAIN'
173 ! include 'COMMON.GEO'
174 ! include 'COMMON.SBRIDGE'
175 ! include 'COMMON.CONTROL'
176 character(len=7) :: chtmp
177 character(len=40) :: chfrm
180 integer :: nft,k,l,i,j,jlee
182 open(icsa_bank,file=csa_bank,status="unknown")
183 write (icsa_bank,900) jlee,nbank,nstep,nft,icycle,cutdif
184 write (icsa_bank,902) nglob_csa, eglob_csa
185 open (igeom,file=intname,status='UNKNOWN')
187 write (icsa_bank,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
188 if (vdisulf) write (icsa_bank,'(101i4)') &
189 bvar_nss(k),((bvar_ss(j,i,k),j=1,2),i=1,bvar_nss(k))
192 write (icsa_bank,850) (rad2deg*bvar(i,l,j,k),i=1,4)
195 if (bvar_nss(k).le.9) then
196 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
197 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
199 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
200 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
201 write (igeom,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
202 bvar_ss(2,i,k),i=10,bvar_nss(k))
204 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
205 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
206 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
207 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
212 if (nstep/200.gt.ilastnstep) then
214 ilastnstep=(ilastnstep+1)*1.5
215 write(chfrm,'(a2,i1,a1)') '(i',int(dlog10(dble(nstep))+1),')'
216 write(chtmp,chfrm) nstep
217 open(icsa_int,file=prefix(:ilen(prefix)) &
218 //'_'//chtmp(:ilen(chtmp))//'.int',status='UNKNOWN')
220 if (bvar_nss(k).le.9) then
221 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
222 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
224 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
225 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
226 write (icsa_int,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
227 bvar_ss(2,i,k),i=10,bvar_nss(k))
229 write (icsa_int,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
230 write (icsa_int,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
231 write (icsa_int,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
232 write (icsa_int,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
240 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
242 902 format (1x,'nglob_csa =',i4,' eglob_csa =',1pe14.5)
243 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
247 end subroutine write_bank
248 !-----------------------------------------------------------------------------
249 subroutine write_bank_reminimized(jlee,nft)
252 use geometry_data, only: nres,rad2deg
254 ! implicit real*8 (a-h,o-z)
255 ! include 'DIMENSIONS'
256 ! include 'COMMON.IOUNITS'
257 ! include 'COMMON.CSA'
258 ! include 'COMMON.BANK'
259 ! include 'COMMON.CHAIN'
260 ! include 'COMMON.GEO'
261 ! include 'COMMON.SBRIDGE'
263 integer :: nft,i,l,j,k,jlee
265 open(icsa_bank_reminimized,file=csa_bank_reminimized,&
267 write (icsa_bank_reminimized,900) &
268 jlee,nbank,nstep,nft,icycle,cutdif
269 open (igeom,file=intname,status='UNKNOWN')
271 write (icsa_bank_reminimized,952) k,bene(k),brmsn(k),&
275 write (icsa_bank_reminimized,850) (rad2deg*bvar(i,l,j,k),i=1,4)
279 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
280 nss,(ihpb(i),jhpb(i),i=1,nss)
282 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
283 nss,(ihpb(i),jhpb(i),i=1,9)
284 write (igeom,'(3X,11(1X,2I3))') (ihpb(i),jhpb(i),i=10,nss)
286 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
287 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
288 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
289 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
291 close(icsa_bank_reminimized)
296 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
298 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
302 end subroutine write_bank_reminimized
303 !-----------------------------------------------------------------------------
304 subroutine read_bank(jlee,nft,cutdifr)
307 use control_data, only: vdisulf
308 use geometry_data, only: nres,deg2rad
310 ! implicit real*8 (a-h,o-z)
311 ! include 'DIMENSIONS'
312 ! include 'COMMON.IOUNITS'
313 ! include 'COMMON.CSA'
314 ! include 'COMMON.BANK'
315 ! include 'COMMON.CHAIN'
316 ! include 'COMMON.GEO'
317 ! include 'COMMON.CONTROL'
318 ! include 'COMMON.SBRIDGE'
319 character(len=80) :: karta
323 integer :: nft,kk,k,l,i,j,jlee
324 real(kind=8) :: cutdifr
326 open(icsa_bank,file=csa_bank,status="old")
327 read (icsa_bank,901) jlee,nbank,nstep,nft,icycle,cutdifr
328 read (icsa_bank,902) nglob_csa, eglob_csa
329 ! if(jleer.ne.jlee) then
330 ! write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,
332 ! call mpi_abort(mpi_comm_world,ierror,ierrcode)
337 read (icsa_bank,'(a80)') karta
338 write(iout,*) "READ_BANK: kk=",kk
340 ! if (index(karta,"*").gt.0) then
341 ! write (iout,*) "***** Stars in bank ***** k=",k,
345 ! read (33,850) (rdummy,i=1,4)
350 call reada(karta,"total E",bene(kk),1.0d20)
351 call reada(karta,"rmsd from N",brmsn(kk),0.0d0)
352 call reada(karta,"%NC",bpncn(kk),0.0d0)
353 read (karta(ilen(karta)-1:),*,end=111,err=111) ibank(kk)
357 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
358 "%NC",bpncn(kk),ibank(kk)
359 ! read (icsa_bank,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
361 read (icsa_bank,'(101i4)') &
362 bvar_nss(kk),((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
363 bvar_ns(kk)=ns-2*bvar_nss(kk)
364 write(iout,*) 'read SSBOND',bvar_nss(kk),&
365 ((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
366 !d write(iout,*) 'read CYS #free ', bvar_ns(kk)
370 do while( iss(i).ne.bvar_ss(1,j,kk)-nres .and. &
371 iss(i).ne.bvar_ss(2,j,kk)-nres .and. &
375 if (j.gt.bvar_nss(kk)) then
380 !d write(iout,*)'read CYS free',(bvar_s(l,kk),l=1,bvar_ns(kk))
384 read (icsa_bank,850) (bvar(i,l,j,kk),i=1,4)
385 ! write (iout,850) (bvar(i,l,j,kk),i=1,4)
387 bvar(i,l,j,kk)=deg2rad*bvar(i,l,j,kk)
394 if (kk.lt.nbank) nbank=kk
395 !d write (*,*) "read_bank ******************* kk",kk,
401 !d write (*,850) (bvar(i,l,j,kk),i=1,4)
407 ! read (33,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
410 ! read (33,850) (bvar(i,l,j,k),i=1,4)
412 ! bvar(i,l,j,k)=deg2rad*bvar(i,l,j,k)
420 952 format (1x,'#',i4,' total E ',f12.3,' rmsd from N ',f8.3,i5)
421 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
422 902 format (1x,11x,i4,12x,1pe14.5)
423 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2,i5)
426 end subroutine read_bank
427 !-----------------------------------------------------------------------------
428 subroutine write_bank1(jlee)
431 use geometry_data, only: nres,rad2deg
432 ! implicit real*8 (a-h,o-z)
433 ! include 'DIMENSIONS'
434 ! include 'COMMON.IOUNITS'
435 ! include 'COMMON.CSA'
436 ! include 'COMMON.BANK'
437 ! include 'COMMON.CHAIN'
438 ! include 'COMMON.GEO'
440 integer :: k,i,l,j,jlee
442 #if defined(AIX) || defined(PGI)
443 open(icsa_bank1,file=csa_bank1,position="append")
445 open(icsa_bank1,file=csa_bank1,access="append")
447 write (icsa_bank1,900) jlee,nbank,nstep,cutdif
449 write (icsa_bank1,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
452 write (icsa_bank1,850) (rad2deg*bvar(i,l,j,k),i=1,4)
458 900 format (4x,"jlee =",i5,3x,"nbank =",i5,3x,"nstep =",i10,f15.5)
459 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
463 end subroutine write_bank1
464 !-----------------------------------------------------------------------------
466 !-----------------------------------------------------------------------------
467 ! subroutine cartprint
469 ! use geometry_data, only: c
470 ! use energy_data, only: itype
471 ! implicit real*8 (a-h,o-z)
472 ! include 'DIMENSIONS'
473 ! include 'COMMON.CHAIN'
474 ! include 'COMMON.INTERACT'
475 ! include 'COMMON.NAMES'
476 ! include 'COMMON.IOUNITS'
481 ! write (iout,110) restyp(itype(i,1)),i,c(1,i),c(2,i),&
482 ! c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
484 ! 100 format (//' alpha-carbon coordinates ',&
485 ! ' centroid coordinates'/ &
486 ! ' ', 6X,'X',11X,'Y',11X,'Z',&
487 ! 10X,'X',11X,'Y',11X,'Z')
488 ! 110 format (a,'(',i3,')',6f12.5)
490 ! end subroutine cartprint
491 !-----------------------------------------------------------------------------
493 !-----------------------------------------------------------------------------
494 subroutine secstrp2dihc
498 ! implicit real*8 (a-h,o-z)
499 ! include 'DIMENSIONS'
500 ! include 'COMMON.GEO'
501 ! include 'COMMON.BOUNDS'
502 ! include 'COMMON.CHAIN'
503 ! include 'COMMON.TORCNSTR'
504 ! include 'COMMON.IOUNITS'
505 !el character(len=1),dimension(nres) :: secstruc !(maxres)
506 !el COMMON/SECONDARYS/secstruc
507 character(len=80) :: line
512 integer :: i,ii,lenpre
514 allocate(secstruc(nres))
516 !dr call getenv_loc('SECPREDFIL',secpred)
518 secpred=prefix(:lenpre)//'.spred'
520 #if defined(WINIFL) || defined(WINPGI)
521 open(isecpred,file=secpred,status='old',readonly,shared)
522 #elif (defined CRAY) || (defined AIX)
523 open(isecpred,file=secpred,status='old',action='read')
525 open(isecpred,file=secpred,status='old')
527 open(isecpred,file=secpred,status='old',action='read')
529 ! read secondary structure prediction from JPRED here!
530 ! read(isecpred,'(A80)',err=100,end=100) line
531 ! read(line,'(f10.3)',err=110) ftors
532 read(isecpred,'(f10.3)',err=110) ftors(1)
534 write (iout,*) 'FTORS factor =',ftors(1)
535 ! initialize secstruc to any
542 call read_secstr_pred(isecpred,iout,errflag)
544 write(iout,*)'There is a problem with the list of secondary-',&
545 'structure prediction'
548 ! 8/13/98 Set limits to generating the dihedral angles
557 if ( secstruc(i) .eq. 'H') then
558 ! Helix restraints for this residue
561 phi0(ii) = 45.0D0*deg2rad
562 drange(ii)= 5.0D0*deg2rad
563 phibound(1,i) = phi0(ii)-drange(ii)
564 phibound(2,i) = phi0(ii)+drange(ii)
565 else if (secstruc(i) .eq. 'E') then
566 ! strand restraints for this residue
569 phi0(ii) = 180.0D0*deg2rad
570 drange(ii)= 5.0D0*deg2rad
571 phibound(1,i) = phi0(ii)-drange(ii)
572 phibound(2,i) = phi0(ii)+drange(ii)
574 ! no restraints for this residue
575 ndih_nconstr=ndih_nconstr+1
576 idih_nconstr(ndih_nconstr)=i
580 ! deallocate(secstruc)
583 write(iout,'(A30,A80)')'Error reading file SECPRED',secpred
584 ! deallocate(secstruc)
587 write(iout,'(A20)')'Error reading FTORS'
588 ! deallocate(secstruc)
590 end subroutine secstrp2dihc
591 !-----------------------------------------------------------------------------
592 subroutine read_secstr_pred(jin,jout,errors)
594 ! implicit real*8 (a-h,o-z)
595 ! INCLUDE 'DIMENSIONS'
596 ! include 'COMMON.IOUNITS'
597 ! include 'COMMON.CHAIN'
598 !el character(len=1),dimension(nres) :: secstruc !(maxres)
599 !el COMMON/SECONDARYS/secstruc
601 character(len=80) :: line,line1 !,ucase
602 logical :: errflag,errors,blankline
605 integer :: jin,jout,iseq,ipos,ipos1,iend,il,&
608 read (jin,'(a)') line
609 write (jout,'(2a)') '> ',line(1:78)
611 ! Remember that we number full residues starting from 2, then, iseq=1 and iseq=nres
612 ! correspond to the end-groups. ADD to the secondary structure prediction "-" for the
613 ! end-groups in the input file "*.spred"
616 do while (index(line1,'$END').eq.0)
617 ! Override commented lines.
620 do while (.not.blankline)
622 call mykey(line,line1,ipos,blankline,errflag)
623 if (errflag) write (jout,'(2a)') &
624 'Error when reading sequence in line: ',line
625 errors=errors .or. errflag
626 if (.not. blankline .and. .not. errflag) then
629 !el if (iseq.le.maxres) then
630 if (line1(1:1).eq.'-' ) then
631 secstruc(iseq)=line1(1:1)
632 else if ( ( ucase(line1(1:1)).eq.'E' ) .or. &
633 ( ucase(line1(1:1)).eq.'H' ) ) then
634 secstruc(iseq)=ucase(line1(1:1))
637 write (jout,1010) line1(1:1), iseq
642 !el write (jout,1000) iseq,maxres
645 do while (ipos1.le.iend)
650 !el if (iseq.le.maxres) then
651 if (line1(ipos1-1:ipos1-1).eq.'-' ) then
652 secstruc(iseq)=line1(ipos1-1:ipos1-1)
653 else if((ucase(line1(ipos1-1:ipos1-1)).eq.'E').or. &
654 (ucase(line1(ipos1-1:ipos1-1)).eq.'H') ) then
655 secstruc(iseq)=ucase(line1(ipos1-1:ipos1-1))
658 write (jout,1010) line1(ipos1-1:ipos1-1), iseq
663 !el write (jout,1000) iseq,maxres
670 read (jin,'(a)') line
671 write (jout,'(2a)') '> ',line(1:78)
675 !d write (jout,'(10a8)') (sequence(i),i=1,iseq-1)
677 !d check whether the found length of the chain is correct.
678 length_of_chain=iseq-1
679 if (length_of_chain .ne. nres) then
681 write (jout,'(a,i4,a,i4,a)') &
682 'Error: the number of labels specified in $SEC_STRUC_PRED (' &
683 ,length_of_chain,') does not match with the number of residues (' &
688 1000 format('Error - the number of residues (',i4,&
689 ') has exceeded maximum (',i4,').')
690 1010 format ('Error - unrecognized secondary structure label',a4,&
693 end subroutine read_secstr_pred
695 !-----------------------------------------------------------------------------
697 !-----------------------------------------------------------------------------
702 use control_data, only:maxterm !,maxtor
706 use control, only: getenv_loc
708 ! Read the parameters of the probability distributions of the virtual-bond
709 ! valence angles and the side chains and energy parameters.
711 ! Important! Energy-term weights ARE NOT read here; they are read from the
712 ! main input file instead, because NO defaults have yet been set for these
715 ! implicit real*8 (a-h,o-z)
716 ! include 'DIMENSIONS'
721 ! include 'COMMON.IOUNITS'
722 ! include 'COMMON.CHAIN'
723 ! include 'COMMON.INTERACT'
724 ! include 'COMMON.GEO'
725 ! include 'COMMON.LOCAL'
726 ! include 'COMMON.TORSION'
727 ! include 'COMMON.SCCOR'
728 ! include 'COMMON.SCROT'
729 ! include 'COMMON.FFIELD'
730 ! include 'COMMON.NAMES'
731 ! include 'COMMON.SBRIDGE'
732 ! include 'COMMON.MD'
733 ! include 'COMMON.SETUP'
734 character(len=1) :: t1,t2,t3
735 character(len=1) :: onelett(4) = (/"G","A","P","D"/)
736 character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
737 logical :: lprint,LaTeX,SPLIT_FOURIERTOR
738 real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob)
739 real(kind=8),dimension(13) :: buse
740 character(len=3) :: lancuch !,ucase
742 integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm,jj
743 integer :: maxinter,junk,kk,ii,ncatprotparm,nkcctyp
744 real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
745 dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
746 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
747 res1,epsijlip,epspeptube,epssctube,sigmapeptube, &
748 sigmasctube,krad2,ract
749 integer :: ichir1,ichir2,ijunk,irdiff
751 character*80 temp1,mychar
752 ! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
753 !el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
754 !el allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
756 ! For printing parameters after they are read set the following in the UNRES
759 ! setenv PRINT_PARM YES
761 ! To print parameters in LaTeX format rather than as ASCII tables:
765 call getenv_loc("PRINT_PARM",lancuch)
766 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
767 call getenv_loc("LATEX",lancuch)
768 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
770 dwa16=2.0d0**(1.0d0/6.0d0)
772 ! Assign virtual-bond length
775 vblinv2=vblinv*vblinv
777 ! Read the virtual-bond parameters, masses, and moments of inertia
778 ! and Stokes' radii of the peptide group and side chains
780 allocate(dsc(ntyp1)) !(ntyp1)
781 allocate(dsc_inv(ntyp1)) !(ntyp1)
782 allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp)
783 allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
784 allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
785 allocate(nbondterm(ntyp)) !(ntyp)
786 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
787 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
788 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
789 allocate(long_r_sidechain(ntyp))
790 allocate(short_r_sidechain(ntyp))
795 allocate(msc(ntyp+1)) !(ntyp+1)
796 allocate(isc(ntyp+1)) !(ntyp+1)
797 allocate(restok(ntyp+1)) !(ntyp+1)
799 read (ibond,*) vbldp0,akp,mp,ip,pstok
802 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
803 dsc(i) = vbldsc0(1,i)
807 dsc_inv(i)=1.0D0/dsc(i)
816 allocate(msc(-ntyp-1:ntyp+1,5)) !(ntyp+1)
817 allocate(isc(-ntyp-1:ntyp+1,5)) !(ntyp+1)
818 allocate(restok(-ntyp-1:ntyp+1,5)) !(ntyp+1)
820 read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp(1),ip(1),pstok(1)
822 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
823 j=1,nbondterm(i)),msc(i,1),isc(i,1),restok(i,1)
824 dsc(i) = vbldsc0(1,i)
828 dsc_inv(i)=1.0D0/dsc(i)
834 read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
837 read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),isc(i,2),restok(i,2)
838 ! dsc(i) = vbldsc0_nucl(1,i)
842 ! dsc_inv(i)=1.0D0/dsc(i)
845 ! read (ibond_nucl,*) junk,vbldp0_nucl,akp_nucl,rjunk,mp(2),ip(2),pstok(2)
846 ! do i=1,ntyp_molec(2)
847 ! read (ibond_nucl,*) nbondterm_nucl(i),(vbldsc0_nucl(j,i),&
848 ! aksc_nucl(j,i),abond0_nucl(j,i),&
849 ! j=1,nbondterm_nucl(i)),msc(i,2),isc(i,2),restok(i,2)
850 ! dsc(i) = vbldsc0(1,i)
854 ! dsc_inv(i)=1.0D0/dsc(i)
859 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
860 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
862 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp(1),ip(1),pstok(1)
864 write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(i),&
865 vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i,1),isc(i,1),restok(i,1)
867 write (iout,'(13x,3f10.5)') &
868 vbldsc0(j,i),aksc(j,i),abond0(j,i)
875 if (.not.allocated(ichargecat)) &
876 allocate (ichargecat(-ntyp_molec(5):ntyp_molec(5)))
878 if (oldion.eq.1) then
880 read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
881 print *,msc(i,5),restok(i,5)
885 read (iion,*) ncatprotparm
886 allocate(catprm(ncatprotparm,4))
888 read (iion,*) (catprm(i,k),i=1,ncatprotparm)
892 allocate(catnuclprm(14,ntyp_molec(2),ntyp_molec(5)))
896 read(iionnucl,*) (catnuclprm(k,j,i),k=1,14)
899 write(*,'(3(5x,a6)11(7x,a6))') "w1 ","w2 ","epslj ","pis1 ", &
900 "sigma0","epsi0 ","chi1 ","chip1 ","sig ","b1 ","b2 ", &
904 write(*,'(3(f10.3,x),11(f12.6,x),a3,2a)') (catnuclprm(k,j,i),k=1,14), &
905 restyp(i,5),"-",restyp(j,2)
908 ! read (iion,*) (vcatprm(k),k=1,ncatprotpram)
909 !----------------------------------------------------
910 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
911 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
912 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
913 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
914 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
915 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
925 allocate(liptranene(ntyp))
926 !C reading lipid parameters
927 write (iout,*) "iliptranpar",iliptranpar
929 read(iliptranpar,*) pepliptran
932 read(iliptranpar,*) liptranene(i)
933 print *,liptranene(i)
936 ! water parmaters entalphy
937 allocate(awaterenta(0:400))
938 allocate(bwaterenta(0:400))
939 allocate(cwaterenta(0:400))
940 allocate(dwaterenta(0:400))
941 allocate(awaterentro(0:400))
942 allocate(bwaterentro(0:400))
943 allocate(cwaterentro(0:400))
944 allocate(dwaterentro(0:400))
946 read(iwaterwater,*) mychar
947 read(iwaterwater,*) ract,awaterenta(0),bwaterenta(0),&
948 cwaterenta(0),dwaterenta(0)
950 read(iwaterwater,*) ract,awaterenta(i),bwaterenta(i),&
951 cwaterenta(i),dwaterenta(i)
952 irdiff=int((ract-2.06d0)*50.0d0)+1
953 if (i.ne.irdiff) print *,"WARTINING",i,ract, irdiff
955 ! water parmaters entrophy
956 read(iwaterwater,*) mychar
957 read(iwaterwater,*) ract,awaterentro(0),bwaterentro(0),&
958 cwaterentro(0),dwaterentro(0)
960 read(iwaterwater,*) ract,awaterentro(i),bwaterentro(i),&
961 cwaterentro(i),dwaterentro(i)
962 irdiff=int((ract-2.06d0)*50.0d0)+1
963 if (i.ne.irdiff) print *,"WARTINING",i,ract, irdiff
969 ! Read the parameters of the probability distribution/energy expression
970 ! of the virtual-bond valence angles theta
973 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
974 (bthet(j,i,1,1),j=1,2)
975 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
976 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
977 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
981 athet(1,i,1,-1)=athet(1,i,1,1)
982 athet(2,i,1,-1)=athet(2,i,1,1)
983 bthet(1,i,1,-1)=-bthet(1,i,1,1)
984 bthet(2,i,1,-1)=-bthet(2,i,1,1)
985 athet(1,i,-1,1)=-athet(1,i,1,1)
986 athet(2,i,-1,1)=-athet(2,i,1,1)
987 bthet(1,i,-1,1)=bthet(1,i,1,1)
988 bthet(2,i,-1,1)=bthet(2,i,1,1)
992 athet(1,i,-1,-1)=athet(1,-i,1,1)
993 athet(2,i,-1,-1)=-athet(2,-i,1,1)
994 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
995 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
996 athet(1,i,-1,1)=athet(1,-i,1,1)
997 athet(2,i,-1,1)=-athet(2,-i,1,1)
998 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
999 bthet(2,i,-1,1)=bthet(2,-i,1,1)
1000 athet(1,i,1,-1)=-athet(1,-i,1,1)
1001 athet(2,i,1,-1)=athet(2,-i,1,1)
1002 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1003 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1004 theta0(i)=theta0(-i)
1008 polthet(j,i)=polthet(j,-i)
1011 gthet(j,i)=gthet(j,-i)
1017 if (.not.LaTeX) then
1018 write (iout,'(a)') &
1019 'Parameters of the virtual-bond valence angles:'
1020 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
1021 ' ATHETA0 ',' A1 ',' A2 ',&
1024 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
1025 a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
1027 write (iout,'(/a/9x,5a/79(1h-))') &
1028 'Parameters of the expression for sigma(theta_c):',&
1029 ' ALPH0 ',' ALPH1 ',' ALPH2 ',&
1030 ' ALPH3 ',' SIGMA0C '
1032 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
1033 (polthet(j,i),j=0,3),sigc0(i)
1035 write (iout,'(/a/9x,5a/79(1h-))') &
1036 'Parameters of the second gaussian:',&
1037 ' THETA0 ',' SIGMA0 ',' G1 ',&
1040 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,theta0(i),&
1041 sig0(i),(gthet(j,i),j=1,3)
1044 write (iout,'(a)') &
1045 'Parameters of the virtual-bond valence angles:'
1046 write (iout,'(/a/9x,5a/79(1h-))') &
1047 'Coefficients of expansion',&
1048 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
1049 ' b1*10^1 ',' b2*10^1 '
1051 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),&
1052 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
1053 (10*bthet(j,i,1,1),j=1,2)
1055 write (iout,'(/a/9x,5a/79(1h-))') &
1056 'Parameters of the expression for sigma(theta_c):',&
1057 ' alpha0 ',' alph1 ',' alph2 ',&
1058 ' alhp3 ',' sigma0c '
1060 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i,1),&
1061 (polthet(j,i),j=0,3),sigc0(i)
1063 write (iout,'(/a/9x,5a/79(1h-))') &
1064 'Parameters of the second gaussian:',&
1065 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
1068 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),theta0(i),&
1069 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
1075 ! Read the parameters of Utheta determined from ab initio surfaces
1076 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
1078 IF (tor_mode.eq.0) THEN
1079 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
1080 ntheterm3,nsingle,ndouble
1081 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
1083 !----------------------------------------------------
1084 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1085 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
1086 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1087 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1088 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
1089 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1090 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1091 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1092 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1093 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1094 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1095 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1096 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1097 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1098 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1099 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1100 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1101 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1102 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1103 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1104 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1105 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1106 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1107 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1109 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
1111 ithetyp(i)=-ithetyp(-i)
1114 aa0thet(:,:,:,:)=0.0d0
1115 aathet(:,:,:,:,:)=0.0d0
1116 bbthet(:,:,:,:,:,:)=0.0d0
1117 ccthet(:,:,:,:,:,:)=0.0d0
1118 ddthet(:,:,:,:,:,:)=0.0d0
1119 eethet(:,:,:,:,:,:)=0.0d0
1120 ffthet(:,:,:,:,:,:,:)=0.0d0
1121 ggthet(:,:,:,:,:,:,:)=0.0d0
1123 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1125 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1126 ! VAR:1=non-glicyne non-proline 2=proline
1127 ! VAR:negative values for D-aminoacid
1129 do j=-nthetyp,nthetyp
1130 do k=-nthetyp,nthetyp
1131 read (ithep,'(6a)',end=111,err=111) res1
1132 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1133 ! VAR: aa0thet is variable describing the average value of Foureir
1134 ! VAR: expansion series
1135 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1136 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1137 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1138 read (ithep,*,end=111,err=111) &
1139 (aathet(l,i,j,k,iblock),l=1,ntheterm)
1140 read (ithep,*,end=111,err=111) &
1141 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1142 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1143 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1144 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1146 read (ithep,*,end=111,err=111) &
1147 (((ffthet(llll,lll,ll,i,j,k,iblock),&
1148 ffthet(lll,llll,ll,i,j,k,iblock),&
1149 ggthet(llll,lll,ll,i,j,k,iblock),&
1150 ggthet(lll,llll,ll,i,j,k,iblock),&
1151 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1156 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1157 ! coefficients of theta-and-gamma-dependent terms are zero.
1158 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1159 ! RECOMENTDED AFTER VERSION 3.3)
1163 ! aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1164 ! aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1166 ! aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1167 ! aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1170 ! aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1172 ! aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1175 ! AND COMMENT THE LOOPS BELOW
1179 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1180 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1182 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1183 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1186 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1188 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1193 ! Substitution for D aminoacids from symmetry.
1196 do j=-nthetyp,nthetyp
1197 do k=-nthetyp,nthetyp
1198 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1200 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1204 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1205 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1206 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1207 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1213 ffthet(llll,lll,ll,i,j,k,iblock)= &
1214 ffthet(llll,lll,ll,-i,-j,-k,iblock)
1215 ffthet(lll,llll,ll,i,j,k,iblock)= &
1216 ffthet(lll,llll,ll,-i,-j,-k,iblock)
1217 ggthet(llll,lll,ll,i,j,k,iblock)= &
1218 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1219 ggthet(lll,llll,ll,i,j,k,iblock)= &
1220 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1229 ! Control printout of the coefficients of virtual-bond-angle potentials
1232 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1237 write (iout,'(//4a)') &
1238 'Type ',onelett(i),onelett(j),onelett(k)
1239 write (iout,'(//a,10x,a)') " l","a[l]"
1240 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1241 write (iout,'(i2,1pe15.5)') &
1242 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1244 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1245 "b",l,"c",l,"d",l,"e",l
1247 write (iout,'(i2,4(1pe15.5))') m,&
1248 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1249 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1253 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1254 "f+",l,"f-",l,"g+",l,"g-",l
1257 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1258 ffthet(n,m,l,i,j,k,iblock),&
1259 ffthet(m,n,l,i,j,k,iblock),&
1260 ggthet(n,m,l,i,j,k,iblock),&
1261 ggthet(m,n,l,i,j,k,iblock)
1272 !C here will be the apropriate recalibrating for D-aminoacid
1273 read (ithep,*,end=121,err=121) nthetyp
1274 allocate(nbend_kcc_Tb(-nthetyp:nthetyp))
1275 allocate(v1bend_chyb(0:36,-nthetyp:nthetyp))
1276 do i=-nthetyp+1,nthetyp-1
1277 read (ithep,*,end=121,err=121) nbend_kcc_Tb(i)
1278 do j=0,nbend_kcc_Tb(i)
1279 read (ithep,*,end=121,err=121) ijunk,v1bend_chyb(j,i)
1283 write (iout,'(a)') &
1284 "Parameters of the valence-only potentials"
1285 do i=-nthetyp+1,nthetyp-1
1286 write (iout,'(2a)') "Type ",toronelet(i)
1287 do k=0,nbend_kcc_Tb(i)
1288 write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
1294 write (2,*) "Start reading THETA_PDB",ithep_pdb
1296 ! write (2,*) 'i=',i
1297 read (ithep_pdb,*,err=111,end=111) &
1298 a0thet(i),(athet(j,i,1,1),j=1,2),&
1299 (bthet(j,i,1,1),j=1,2)
1300 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1301 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1302 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1303 sigc0(i)=sigc0(i)**2
1306 athet(1,i,1,-1)=athet(1,i,1,1)
1307 athet(2,i,1,-1)=athet(2,i,1,1)
1308 bthet(1,i,1,-1)=-bthet(1,i,1,1)
1309 bthet(2,i,1,-1)=-bthet(2,i,1,1)
1310 athet(1,i,-1,1)=-athet(1,i,1,1)
1311 athet(2,i,-1,1)=-athet(2,i,1,1)
1312 bthet(1,i,-1,1)=bthet(1,i,1,1)
1313 bthet(2,i,-1,1)=bthet(2,i,1,1)
1316 a0thet(i)=a0thet(-i)
1317 athet(1,i,-1,-1)=athet(1,-i,1,1)
1318 athet(2,i,-1,-1)=-athet(2,-i,1,1)
1319 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
1320 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
1321 athet(1,i,-1,1)=athet(1,-i,1,1)
1322 athet(2,i,-1,1)=-athet(2,-i,1,1)
1323 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
1324 bthet(2,i,-1,1)=bthet(2,-i,1,1)
1325 athet(1,i,1,-1)=-athet(1,-i,1,1)
1326 athet(2,i,1,-1)=athet(2,-i,1,1)
1327 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1328 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1329 theta0(i)=theta0(-i)
1333 polthet(j,i)=polthet(j,-i)
1336 gthet(j,i)=gthet(j,-i)
1339 write (2,*) "End reading THETA_PDB"
1343 !--------------- Reading theta parameters for nucleic acid-------
1344 read (ithep_nucl,*,err=111,end=111) nthetyp_nucl,ntheterm_nucl,&
1345 ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
1346 nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
1347 allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1348 allocate(aa0thet_nucl(nthetyp_nucl+1,&
1349 nthetyp_nucl+1,nthetyp_nucl+1))
1350 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1351 allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
1352 nthetyp_nucl+1,nthetyp_nucl+1))
1353 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1354 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1355 allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1356 nthetyp_nucl+1,nthetyp_nucl+1))
1357 allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1358 nthetyp_nucl+1,nthetyp_nucl+1))
1359 allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1360 nthetyp_nucl+1,nthetyp_nucl+1))
1361 allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1362 nthetyp_nucl+1,nthetyp_nucl+1))
1363 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1364 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1365 allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1366 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1367 allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1368 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1370 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1371 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1373 read (ithep_nucl,*,err=111,end=111) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
1375 aa0thet_nucl(:,:,:)=0.0d0
1376 aathet_nucl(:,:,:,:)=0.0d0
1377 bbthet_nucl(:,:,:,:,:)=0.0d0
1378 ccthet_nucl(:,:,:,:,:)=0.0d0
1379 ddthet_nucl(:,:,:,:,:)=0.0d0
1380 eethet_nucl(:,:,:,:,:)=0.0d0
1381 ffthet_nucl(:,:,:,:,:,:)=0.0d0
1382 ggthet_nucl(:,:,:,:,:,:)=0.0d0
1387 read (ithep_nucl,'(3a)',end=111,err=111) t1,t2,t3
1388 read (ithep_nucl,*,end=111,err=111) aa0thet_nucl(i,j,k)
1389 read (ithep_nucl,*,end=111,err=111)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
1390 read (ithep_nucl,*,end=111,err=111) &
1391 (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1392 (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1393 (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1394 (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
1395 read (ithep_nucl,*,end=111,err=111) &
1396 (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
1397 ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
1398 llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
1403 !-------------------------------------------
1404 allocate(nlob(ntyp1)) !(ntyp1)
1405 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1406 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1407 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1414 gaussc(:,:,:,:)=0.0D0
1418 ! Read the parameters of the probability distribution/energy expression
1419 ! of the side chains.
1422 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1426 dsc_inv(i)=1.0D0/dsc(i)
1437 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1438 ((blower(k,l,1),l=1,k),k=1,3)
1439 censc(1,1,-i)=censc(1,1,i)
1440 censc(2,1,-i)=censc(2,1,i)
1441 censc(3,1,-i)=-censc(3,1,i)
1443 read (irotam,*,end=112,err=112) bsc(j,i)
1444 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1445 ((blower(k,l,j),l=1,k),k=1,3)
1446 censc(1,j,-i)=censc(1,j,i)
1447 censc(2,j,-i)=censc(2,j,i)
1448 censc(3,j,-i)=-censc(3,j,i)
1449 ! BSC is amplitude of Gaussian
1456 akl=akl+blower(k,m,j)*blower(l,m,j)
1460 if (((k.eq.3).and.(l.ne.3)) &
1461 .or.((l.eq.3).and.(k.ne.3))) then
1462 gaussc(k,l,j,-i)=-akl
1463 gaussc(l,k,j,-i)=-akl
1465 gaussc(k,l,j,-i)=akl
1466 gaussc(l,k,j,-i)=akl
1475 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1478 if (nlobi.gt.0) then
1480 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i,1),&
1481 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1482 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1483 'log h',(bsc(j,i),j=1,nlobi)
1484 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1485 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1487 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1488 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1491 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1492 write (iout,'(a,f10.4,4(16x,f10.4))') &
1493 'Center ',(bsc(j,i),j=1,nlobi)
1494 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1503 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1504 ! added by Urszula Kozlowska 07/11/2007
1506 !el Maximum number of SC local term fitting function coefficiants
1507 !el integer,parameter :: maxsccoef=65
1509 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1512 read (irotam,*,end=112,err=112)
1514 read (irotam,*,end=112,err=112)
1517 read(irotam,*,end=112,err=112) sc_parmin(j,i)
1522 allocate(nterm_scend(2,ntyp))
1523 allocate(arotam_end(0:6,2,ntyp))
1526 read (irotam_end,*) ijunk
1527 !c write (iout,*) "ijunk",ijunk
1531 read (irotam_end,'(a)')
1532 read (irotam_end,*) nterm_scend(j,i)
1533 !c write (iout,*) "i",i," j",j," nterm",nterm_scend(j,i)
1534 do k=0,nterm_scend(j,i)
1535 read (irotam_end,*) ijunk,arotam_end(k,j,i)
1536 !c write (iout,*) "k",k," arotam",arotam_end(k,j,i)
1542 write (iout,'(a)') &
1543 "Parameters of the local potentials of sidechain ends"
1545 write (iout,'(5x,9x,2hp-,a3,6x,9x,a3,2h-p)')&
1546 restyp(i,1),restyp(i,1)
1547 do j=0,max0(nterm_scend(1,i),nterm_scend(2,i))
1548 write (iout,'(i5,2f20.10)') &
1549 j,arotam_end(j,1,i),arotam_end(j,2,i)
1556 !---------reading nucleic acid parameters for rotamers-------------------
1557 allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp)
1558 do i=1,ntyp_molec(2)
1559 read (irotam_nucl,*,end=112,err=112)
1561 read(irotam_nucl,*,end=112,err=112) sc_parmin_nucl(j,i)
1567 write (iout,*) "Base rotamer parameters"
1568 do i=1,ntyp_molec(2)
1569 write (iout,'(a)') restyp(i,2)
1570 write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
1575 ! Read the parameters of the probability distribution/energy expression
1576 ! of the side chains.
1578 write (2,*) "Start reading ROTAM_PDB"
1580 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1584 dsc_inv(i)=1.0D0/dsc(i)
1595 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1596 ((blower(k,l,1),l=1,k),k=1,3)
1598 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1599 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1600 ((blower(k,l,j),l=1,k),k=1,3)
1607 akl=akl+blower(k,m,j)*blower(l,m,j)
1617 write (2,*) "End reading ROTAM_PDB"
1623 !C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1624 !C interaction energy of the Gly, Ala, and Pro prototypes.
1626 read (ifourier,*) nloctyp
1627 SPLIT_FOURIERTOR = nloctyp.lt.0
1628 nloctyp = iabs(nloctyp)
1629 !C allocate(b1(2,nres)) !(2,-maxtor:maxtor)
1630 !C allocate(b2(2,nres)) !(2,-maxtor:maxtor)
1631 !C allocate(b1tilde(2,nres)) !(2,-maxtor:maxtor)
1632 !C allocate(ctilde(2,2,nres))
1633 !C allocate(dtilde(2,2,nres)) !(2,2,-maxtor:maxtor)
1634 !C allocate(gtb1(2,nres))
1635 !C allocate(gtb2(2,nres))
1636 !C allocate(cc(2,2,nres))
1637 !C allocate(dd(2,2,nres))
1638 !C allocate(ee(2,2,nres))
1639 !C allocate(gtcc(2,2,nres))
1640 !C allocate(gtdd(2,2,nres))
1641 !C allocate(gtee(2,2,nres))
1644 allocate(itype2loc(-ntyp1:ntyp1))
1645 allocate(iloctyp(-nloctyp:nloctyp))
1646 allocate(bnew1(3,2,-nloctyp:nloctyp))
1647 allocate(bnew2(3,2,-nloctyp:nloctyp))
1648 allocate(ccnew(3,2,-nloctyp:nloctyp))
1649 allocate(ddnew(3,2,-nloctyp:nloctyp))
1650 allocate(e0new(3,-nloctyp:nloctyp))
1651 allocate(eenew(2,2,2,-nloctyp:nloctyp))
1652 allocate(bnew1tor(3,2,-nloctyp:nloctyp))
1653 allocate(bnew2tor(3,2,-nloctyp:nloctyp))
1654 allocate(ccnewtor(3,2,-nloctyp:nloctyp))
1655 allocate(ddnewtor(3,2,-nloctyp:nloctyp))
1656 allocate(e0newtor(3,-nloctyp:nloctyp))
1657 allocate(eenewtor(2,2,2,-nloctyp:nloctyp))
1670 read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
1671 read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
1672 itype2loc(ntyp1)=nloctyp
1673 iloctyp(nloctyp)=ntyp1
1675 itype2loc(-i)=-itype2loc(i)
1678 allocate(iloctyp(-nloctyp:nloctyp))
1679 allocate(itype2loc(-ntyp1:ntyp1))
1686 iloctyp(-i)=-iloctyp(i)
1688 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1689 !c write (iout,*) "nloctyp",nloctyp,
1690 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1691 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1692 !c write (iout,*) "nloctyp",nloctyp,
1693 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1696 !c write (iout,*) "NEWCORR",i
1697 read (ifourier,*,end=115,err=115)
1700 read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
1703 !c write (iout,*) "NEWCORR BNEW1"
1704 !c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
1707 read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
1710 !c write (iout,*) "NEWCORR BNEW2"
1711 !c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
1713 read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
1714 read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
1716 !c write (iout,*) "NEWCORR CCNEW"
1717 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1719 read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
1720 read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
1722 !c write (iout,*) "NEWCORR DDNEW"
1723 !c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
1727 read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
1731 !c write (iout,*) "NEWCORR EENEW1"
1732 !c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1734 read (ifourier,*,end=115,err=115) e0new(ii,i)
1736 !c write (iout,*) (e0new(ii,i),ii=1,3)
1738 !c write (iout,*) "NEWCORR EENEW"
1741 ccnew(ii,1,i)=ccnew(ii,1,i)/2
1742 ccnew(ii,2,i)=ccnew(ii,2,i)/2
1743 ddnew(ii,1,i)=ddnew(ii,1,i)/2
1744 ddnew(ii,2,i)=ddnew(ii,2,i)/2
1749 bnew1(ii,1,-i)= bnew1(ii,1,i)
1750 bnew1(ii,2,-i)=-bnew1(ii,2,i)
1751 bnew2(ii,1,-i)= bnew2(ii,1,i)
1752 bnew2(ii,2,-i)=-bnew2(ii,2,i)
1755 !c ccnew(ii,1,i)=ccnew(ii,1,i)/2
1756 !c ccnew(ii,2,i)=ccnew(ii,2,i)/2
1757 !c ddnew(ii,1,i)=ddnew(ii,1,i)/2
1758 !c ddnew(ii,2,i)=ddnew(ii,2,i)/2
1759 ccnew(ii,1,-i)=ccnew(ii,1,i)
1760 ccnew(ii,2,-i)=-ccnew(ii,2,i)
1761 ddnew(ii,1,-i)=ddnew(ii,1,i)
1762 ddnew(ii,2,-i)=-ddnew(ii,2,i)
1764 e0new(1,-i)= e0new(1,i)
1765 e0new(2,-i)=-e0new(2,i)
1766 e0new(3,-i)=-e0new(3,i)
1768 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
1769 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
1770 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
1771 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
1775 write (iout,'(a)') "Coefficients of the multibody terms"
1776 do i=-nloctyp+1,nloctyp-1
1777 write (iout,*) "Type: ",onelet(iloctyp(i))
1778 write (iout,*) "Coefficients of the expansion of B1"
1780 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
1782 write (iout,*) "Coefficients of the expansion of B2"
1784 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
1786 write (iout,*) "Coefficients of the expansion of C"
1787 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
1788 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
1789 write (iout,*) "Coefficients of the expansion of D"
1790 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
1791 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
1792 write (iout,*) "Coefficients of the expansion of E"
1793 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
1796 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
1801 IF (SPLIT_FOURIERTOR) THEN
1803 !c write (iout,*) "NEWCORR TOR",i
1804 read (ifourier,*,end=115,err=115)
1807 read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i)
1810 !c write (iout,*) "NEWCORR BNEW1 TOR"
1811 !c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
1814 read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i)
1817 !c write (iout,*) "NEWCORR BNEW2 TOR"
1818 !c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
1820 read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i)
1821 read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i)
1823 !c write (iout,*) "NEWCORR CCNEW TOR"
1824 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1826 read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i)
1827 read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i)
1829 !c write (iout,*) "NEWCORR DDNEW TOR"
1830 !c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
1834 read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i)
1838 !c write (iout,*) "NEWCORR EENEW1 TOR"
1839 !c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1841 read (ifourier,*,end=115,err=115) e0newtor(ii,i)
1843 !c write (iout,*) (e0newtor(ii,i),ii=1,3)
1845 !c write (iout,*) "NEWCORR EENEW TOR"
1848 ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
1849 ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
1850 ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
1851 ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
1856 bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
1857 bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
1858 bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
1859 bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
1862 ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
1863 ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
1864 ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
1865 ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
1867 e0newtor(1,-i)= e0newtor(1,i)
1868 e0newtor(2,-i)=-e0newtor(2,i)
1869 e0newtor(3,-i)=-e0newtor(3,i)
1871 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
1872 eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
1873 eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
1874 eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
1878 write (iout,'(a)') &
1879 "Single-body coefficients of the torsional potentials"
1880 do i=-nloctyp+1,nloctyp-1
1881 write (iout,*) "Type: ",onelet(iloctyp(i))
1882 write (iout,*) "Coefficients of the expansion of B1tor"
1884 write (iout,'(3hB1(,i1,1h),3f10.5)') &
1885 j,(bnew1tor(k,j,i),k=1,3)
1887 write (iout,*) "Coefficients of the expansion of B2tor"
1889 write (iout,'(3hB2(,i1,1h),3f10.5)') &
1890 j,(bnew2tor(k,j,i),k=1,3)
1892 write (iout,*) "Coefficients of the expansion of Ctor"
1893 write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
1894 write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
1895 write (iout,*) "Coefficients of the expansion of Dtor"
1896 write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
1897 write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
1898 write (iout,*) "Coefficients of the expansion of Etor"
1899 write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
1902 write (iout,'(1hE,2i1,2f10.5)') &
1903 j,k,(eenewtor(l,j,k,i),l=1,2)
1909 do i=-nloctyp+1,nloctyp-1
1912 bnew1tor(ii,j,i)=bnew1(ii,j,i)
1917 bnew2tor(ii,j,i)=bnew2(ii,j,i)
1921 ccnewtor(ii,1,i)=ccnew(ii,1,i)
1922 ccnewtor(ii,2,i)=ccnew(ii,2,i)
1923 ddnewtor(ii,1,i)=ddnew(ii,1,i)
1924 ddnewtor(ii,2,i)=ddnew(ii,2,i)
1927 ENDIF !SPLIT_FOURIER_TOR
1929 allocate(ccold(2,2,-nloctyp-1:nloctyp+1))
1930 allocate(ddold(2,2,-nloctyp-1:nloctyp+1))
1931 allocate(eeold(2,2,-nloctyp-1:nloctyp+1))
1932 allocate(b(13,-nloctyp-1:nloctyp+1))
1934 write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
1936 read (ifourier,*,end=115,err=115)
1937 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
1939 write (iout,*) 'Type ',onelet(iloctyp(i))
1940 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
1954 !c B1tilde(1,i) = b(3)
1955 !c! B1tilde(2,i) =-b(5)
1956 !c! B1tilde(1,-i) =-b(3)
1957 !c! B1tilde(2,-i) =b(5)
1958 !c! b1tilde(1,i)=0.0d0
1959 !c b1tilde(2,i)=0.0d0
1964 !cc B1tilde(1,i) = b(3,i)
1965 !cc B1tilde(2,i) =-b(5,i)
1966 !c B1tilde(1,-i) =-b(3,i)
1967 !c B1tilde(2,-i) =b(5,i)
1968 !cc b1tilde(1,i)=0.0d0
1969 !cc b1tilde(2,i)=0.0d0
1970 !cc B2(1,i) = b(2,i)
1971 !cc B2(2,i) = b(4,i)
1973 !c B2(2,-i) =-b(4,i)
1977 CCold(1,1,i)= b(7,i)
1978 CCold(2,2,i)=-b(7,i)
1979 CCold(2,1,i)= b(9,i)
1980 CCold(1,2,i)= b(9,i)
1981 CCold(1,1,-i)= b(7,i)
1982 CCold(2,2,-i)=-b(7,i)
1983 CCold(2,1,-i)=-b(9,i)
1984 CCold(1,2,-i)=-b(9,i)
1989 !c Ctilde(1,1,i)= CCold(1,1,i)
1990 !c Ctilde(1,2,i)= CCold(1,2,i)
1991 !c Ctilde(2,1,i)=-CCold(2,1,i)
1992 !c Ctilde(2,2,i)=-CCold(2,2,i)
1997 !c Ctilde(1,1,i)= CCold(1,1,i)
1998 !c Ctilde(1,2,i)= CCold(1,2,i)
1999 !c Ctilde(2,1,i)=-CCold(2,1,i)
2000 !c Ctilde(2,2,i)=-CCold(2,2,i)
2002 !c Ctilde(1,1,i)=0.0d0
2003 !c Ctilde(1,2,i)=0.0d0
2004 !c Ctilde(2,1,i)=0.0d0
2005 !c Ctilde(2,2,i)=0.0d0
2006 DDold(1,1,i)= b(6,i)
2007 DDold(2,2,i)=-b(6,i)
2008 DDold(2,1,i)= b(8,i)
2009 DDold(1,2,i)= b(8,i)
2010 DDold(1,1,-i)= b(6,i)
2011 DDold(2,2,-i)=-b(6,i)
2012 DDold(2,1,-i)=-b(8,i)
2013 DDold(1,2,-i)=-b(8,i)
2018 !c Dtilde(1,1,i)= DD(1,1,i)
2019 !c Dtilde(1,2,i)= DD(1,2,i)
2020 !c Dtilde(2,1,i)=-DD(2,1,i)
2021 !c Dtilde(2,2,i)=-DD(2,2,i)
2023 !c Dtilde(1,1,i)=0.0d0
2024 !c Dtilde(1,2,i)=0.0d0
2025 !c Dtilde(2,1,i)=0.0d0
2026 !c Dtilde(2,2,i)=0.0d0
2027 EEold(1,1,i)= b(10,i)+b(11,i)
2028 EEold(2,2,i)=-b(10,i)+b(11,i)
2029 EEold(2,1,i)= b(12,i)-b(13,i)
2030 EEold(1,2,i)= b(12,i)+b(13,i)
2031 EEold(1,1,-i)= b(10,i)+b(11,i)
2032 EEold(2,2,-i)=-b(10,i)+b(11,i)
2033 EEold(2,1,-i)=-b(12,i)+b(13,i)
2034 EEold(1,2,-i)=-b(12,i)-b(13,i)
2035 write(iout,*) "TU DOCHODZE"
2041 !c ee(2,1,i)=ee(1,2,i)
2046 "Coefficients of the cumulants (independent of valence angles)"
2047 do i=-nloctyp+1,nloctyp-1
2048 write (iout,*) 'Type ',onelet(iloctyp(i))
2050 write(iout,'(2f10.5)') B(3,i),B(5,i)
2052 write(iout,'(2f10.5)') B(2,i),B(4,i)
2055 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
2059 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
2063 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
2072 ! Read torsional parameters in old format
2074 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
2076 read (itorp,*,end=113,err=113) ntortyp,nterm_old
2077 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
2078 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
2080 !el from energy module--------
2081 allocate(v1(nterm_old,ntortyp,ntortyp))
2082 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
2083 !el---------------------------
2088 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
2094 write (iout,'(/a/)') 'Torsional constants:'
2097 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
2098 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
2104 ! Read torsional parameters
2106 IF (TOR_MODE.eq.0) THEN
2107 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
2108 read (itorp,*,end=113,err=113) ntortyp
2109 !el from energy module---------
2110 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2111 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2113 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2114 allocate(vlor2(maxlor,ntortyp,ntortyp))
2115 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
2116 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2118 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
2119 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2120 !el---------------------------
2123 !el---------------------------
2125 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
2127 itortyp(i)=-itortyp(-i)
2129 itortyp(ntyp1)=ntortyp
2130 itortyp(-ntyp1)=-ntortyp
2132 write (iout,*) 'ntortyp',ntortyp
2134 do j=-ntortyp+1,ntortyp-1
2135 read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
2137 nterm(-i,-j,iblock)=nterm(i,j,iblock)
2138 nlor(-i,-j,iblock)=nlor(i,j,iblock)
2141 do k=1,nterm(i,j,iblock)
2142 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
2144 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
2145 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
2146 v0ij=v0ij+si*v1(k,i,j,iblock)
2148 ! write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
2149 ! write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
2150 ! v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
2152 do k=1,nlor(i,j,iblock)
2153 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
2154 vlor2(k,i,j),vlor3(k,i,j)
2155 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
2158 v0(-i,-j,iblock)=v0ij
2164 write (iout,'(/a/)') 'Torsional constants:'
2166 do i=-ntortyp,ntortyp
2167 do j=-ntortyp,ntortyp
2168 write (iout,*) 'ityp',i,' jtyp',j
2169 write (iout,*) 'Fourier constants'
2170 do k=1,nterm(i,j,iblock)
2171 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
2174 write (iout,*) 'Lorenz constants'
2175 do k=1,nlor(i,j,iblock)
2176 write (iout,'(3(1pe15.5))') &
2177 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
2183 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
2185 ! 6/23/01 Read parameters for double torsionals
2187 !el from energy module------------
2188 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2189 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2190 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2191 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2192 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2193 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2194 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2195 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2196 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2197 !---------------------------------
2201 do j=-ntortyp+1,ntortyp-1
2202 do k=-ntortyp+1,ntortyp-1
2203 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
2204 ! write (iout,*) "OK onelett",
2207 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
2208 .or. t3.ne.toronelet(k)) then
2209 write (iout,*) "Error in double torsional parameter file",&
2212 call MPI_Finalize(Ierror)
2214 stop "Error in double torsional parameter file"
2216 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
2217 ntermd_2(i,j,k,iblock)
2218 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
2219 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
2220 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
2221 ntermd_1(i,j,k,iblock))
2222 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
2223 ntermd_1(i,j,k,iblock))
2224 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
2225 ntermd_1(i,j,k,iblock))
2226 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
2227 ntermd_1(i,j,k,iblock))
2228 ! Martix of D parameters for one dimesional foureir series
2229 do l=1,ntermd_1(i,j,k,iblock)
2230 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
2231 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
2232 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
2233 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
2234 ! write(iout,*) "whcodze" ,
2235 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
2237 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
2238 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
2239 v2s(m,l,i,j,k,iblock),&
2240 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
2241 ! Martix of D parameters for two dimesional fourier series
2242 do l=1,ntermd_2(i,j,k,iblock)
2244 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
2245 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
2246 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
2247 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
2256 write (iout,*) 'Constants for double torsionals'
2259 do j=-ntortyp+1,ntortyp-1
2260 do k=-ntortyp+1,ntortyp-1
2261 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
2262 ' nsingle',ntermd_1(i,j,k,iblock),&
2263 ' ndouble',ntermd_2(i,j,k,iblock)
2265 write (iout,*) 'Single angles:'
2266 do l=1,ntermd_1(i,j,k,iblock)
2267 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
2268 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
2269 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
2270 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
2273 write (iout,*) 'Pairs of angles:'
2274 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
2275 do l=1,ntermd_2(i,j,k,iblock)
2276 write (iout,'(i5,20f10.5)') &
2277 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
2280 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
2281 do l=1,ntermd_2(i,j,k,iblock)
2282 write (iout,'(i5,20f10.5)') &
2283 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
2284 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
2294 itype2loc(i)=itortyp(i)
2298 ELSE IF (TOR_MODE.eq.1) THEN
2300 !C read valence-torsional parameters
2301 read (itorp,*,end=121,err=121) ntortyp
2303 write (iout,*) "Valence-torsional parameters read in ntortyp",&
2305 read (itorp,*,end=121,err=121) (itortyp(i),i=1,ntyp)
2306 write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
2309 itype2loc(i)=itortyp(i)
2313 itortyp(i)=-itortyp(-i)
2315 do i=-ntortyp+1,ntortyp-1
2316 do j=-ntortyp+1,ntortyp-1
2317 !C first we read the cos and sin gamma parameters
2318 read (itorp,'(13x,a)',end=121,err=121) string
2319 write (iout,*) i,j,string
2320 read (itorp,*,end=121,err=121) &
2321 nterm_kcc(j,i),nterm_kcc_Tb(j,i)
2322 !C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
2323 do k=1,nterm_kcc(j,i)
2324 do l=1,nterm_kcc_Tb(j,i)
2325 do ll=1,nterm_kcc_Tb(j,i)
2326 read (itorp,*,end=121,err=121) ii,jj,kk, &
2327 v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
2335 !c AL 4/8/16: Calculate coefficients from one-body parameters
2337 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
2338 allocate(nterm_kcc(-ntyp1:ntyp1,-ntyp1:ntyp1))
2339 allocate(nterm_kcc_Tb(-ntyp1:ntyp1,-ntyp1:ntyp1))
2340 allocate(v1_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
2341 allocate(v2_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
2344 print *,i,itortyp(i)
2345 itortyp(i)=itype2loc(i)
2348 "Val-tor parameters calculated from cumulant coefficients ntortyp"&
2350 do i=-ntortyp+1,ntortyp-1
2351 do j=-ntortyp+1,ntortyp-1
2354 do k=1,nterm_kcc_Tb(j,i)
2355 do l=1,nterm_kcc_Tb(j,i)
2356 v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)&
2357 +bnew1tor(k,2,i)*bnew2tor(l,2,j)
2358 v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)&
2359 +bnew1tor(k,2,i)*bnew2tor(l,1,j)
2362 do k=1,nterm_kcc_Tb(j,i)
2363 do l=1,nterm_kcc_Tb(j,i)
2365 v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
2366 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
2367 v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
2368 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
2370 v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
2371 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
2372 v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
2373 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
2377 !c f(theta,gamma)=-(b21(theta)*b11(theta)+b12(theta)*b22(theta))*cos(gamma)-(b22(theta)*b11(theta)+b21(theta)*b12(theta))*sin(gamma)+(c11(theta)*d11(theta)-c12(theta)*d12(theta))*cos(2*gamma)+(c12(theta)*d11(theta)+c11(theta)*d12(theta))*sin(2*gamma)
2381 write (iout,*) "TOR_MODE>1 only with NEWCORR"
2386 if (tor_mode.gt.0 .and. lprint) then
2387 !c Print valence-torsional parameters
2388 write (iout,'(a)') &
2389 "Parameters of the valence-torsional potentials"
2390 do i=-ntortyp+1,ntortyp-1
2391 do j=-ntortyp+1,ntortyp-1
2392 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
2393 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
2394 do k=1,nterm_kcc(j,i)
2395 do l=1,nterm_kcc_Tb(j,i)
2396 do ll=1,nterm_kcc_Tb(j,i)
2397 write (iout,'(3i5,2f15.4)')&
2398 k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
2406 allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
2407 read (itorp_nucl,*,end=113,err=113) ntortyp_nucl
2408 ! print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
2409 !el from energy module---------
2410 allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2411 allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2413 allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2414 allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
2415 allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
2416 allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2418 allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
2419 allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2420 !el---------------------------
2423 !el--------------------
2424 read (itorp_nucl,*,end=113,err=113) &
2425 (itortyp_nucl(i),i=1,ntyp_molec(2))
2426 ! print *,itortyp_nucl(:)
2427 !c write (iout,*) 'ntortyp',ntortyp
2430 read (itorp_nucl,*,end=113,err=113) nterm_nucl(i,j),nlor_nucl(i,j)
2431 ! print *,nterm_nucl(i,j),nlor_nucl(i,j)
2434 do k=1,nterm_nucl(i,j)
2435 read (itorp_nucl,*,end=113,err=113) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
2436 v0ij=v0ij+si*v1_nucl(k,i,j)
2439 do k=1,nlor_nucl(i,j)
2440 read (itorp,*,end=113,err=113) kk,vlor1_nucl(k,i,j),&
2441 vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
2442 v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
2448 ! Read of Side-chain backbone correlation parameters
2449 ! Modified 11 May 2012 by Adasko
2452 read (isccor,*,end=119,err=119) nsccortyp
2454 !el from module energy-------------
2455 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2456 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
2457 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
2458 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
2459 !-----------------------------------
2461 !el from module energy-------------
2462 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
2464 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
2466 isccortyp(i)=-isccortyp(-i)
2468 iscprol=isccortyp(20)
2469 ! write (iout,*) 'ntortyp',ntortyp
2471 !c maxinter is maximum interaction sites
2472 !el from module energy---------
2473 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2474 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
2475 -nsccortyp:nsccortyp))
2476 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
2477 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
2478 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
2479 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
2480 !-----------------------------------
2484 read (isccor,*,end=119,err=119) &
2485 nterm_sccor(i,j),nlor_sccor(i,j)
2491 nterm_sccor(-i,j)=nterm_sccor(i,j)
2492 nterm_sccor(-i,-j)=nterm_sccor(i,j)
2493 nterm_sccor(i,-j)=nterm_sccor(i,j)
2494 do k=1,nterm_sccor(i,j)
2495 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
2497 if (j.eq.iscprol) then
2498 if (i.eq.isccortyp(10)) then
2499 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2500 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2502 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
2503 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
2504 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
2505 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
2506 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2507 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2508 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2509 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2512 if (i.eq.isccortyp(10)) then
2513 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2514 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2516 if (j.eq.isccortyp(10)) then
2517 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
2518 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
2520 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
2521 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2522 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2523 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2524 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2525 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2529 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
2530 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
2531 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
2532 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
2535 do k=1,nlor_sccor(i,j)
2536 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
2537 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2538 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
2539 (1+vlor3sccor(k,i,j)**2)
2541 v0sccor(l,i,j)=v0ijsccor
2542 v0sccor(l,-i,j)=v0ijsccor1
2543 v0sccor(l,i,-j)=v0ijsccor2
2544 v0sccor(l,-i,-j)=v0ijsccor3
2550 !el from module energy-------------
2551 allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
2553 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
2554 ! write (iout,*) 'ntortyp',ntortyp
2556 !c maxinter is maximum interaction sites
2557 !el from module energy---------
2558 allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2559 allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
2560 allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
2561 allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
2562 !-----------------------------------
2566 read (isccor,*,end=119,err=119) &
2567 nterm_sccor(i,j),nlor_sccor(i,j)
2571 do k=1,nterm_sccor(i,j)
2572 print *,"test",i,j,k,l
2573 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
2575 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
2578 do k=1,nlor_sccor(i,j)
2579 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
2580 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2581 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
2582 (1+vlor3sccor(k,i,j)**2)
2584 v0sccor(l,i,j)=v0ijsccor !el ,iblock
2593 write (iout,'(/a/)') 'Torsional constants:'
2596 write (iout,*) 'ityp',i,' jtyp',j
2597 write (iout,*) 'Fourier constants'
2598 do k=1,nterm_sccor(i,j)
2599 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
2601 write (iout,*) 'Lorenz constants'
2602 do k=1,nlor_sccor(i,j)
2603 write (iout,'(3(1pe15.5))') &
2604 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2611 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
2612 ! interaction energy of the Gly, Ala, and Pro prototypes.
2615 ! Read electrostatic-interaction parameters
2620 write (iout,'(/a)') 'Electrostatic interaction constants:'
2621 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
2622 'IT','JT','APP','BPP','AEL6','AEL3'
2624 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
2625 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
2626 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
2627 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
2632 app (i,j)=epp(i,j)*rri*rri
2633 bpp (i,j)=-2.0D0*epp(i,j)*rri
2634 ael6(i,j)=elpp6(i,j)*4.2D0**6
2635 ael3(i,j)=elpp3(i,j)*4.2D0**3
2637 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
2643 ! Read side-chain interaction parameters.
2645 !el from module energy - COMMON.INTERACT-------
2646 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
2647 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
2648 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
2650 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
2651 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
2652 allocate(epslip(ntyp,ntyp))
2660 !--------------------------------
2662 read (isidep,*,end=117,err=117) ipot,expon
2663 if (ipot.lt.1 .or. ipot.gt.5) then
2664 write (iout,'(2a)') 'Error while reading SC interaction',&
2665 'potential file - unknown potential type.'
2667 call MPI_Finalize(Ierror)
2673 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
2674 ', exponents are ',expon,2*expon
2675 ! goto (10,20,30,30,40) ipot
2677 !----------------------- LJ potential ---------------------------------
2679 ! 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2680 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2681 (sigma0(i),i=1,ntyp)
2683 write (iout,'(/a/)') 'Parameters of the LJ potential:'
2684 write (iout,'(a/)') 'The epsilon array:'
2685 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2686 write (iout,'(/a)') 'One-body parameters:'
2687 write (iout,'(a,4x,a)') 'residue','sigma'
2688 write (iout,'(a3,6x,f10.5)') (restyp(i,1),sigma0(i),i=1,ntyp)
2691 !----------------------- LJK potential --------------------------------
2693 ! 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2694 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2695 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
2697 write (iout,'(/a/)') 'Parameters of the LJK potential:'
2698 write (iout,'(a/)') 'The epsilon array:'
2699 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2700 write (iout,'(/a)') 'One-body parameters:'
2701 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
2702 write (iout,'(a3,6x,2f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2706 !---------------------- GB or BP potential -----------------------------
2709 ! print *,"I AM in SCELE",scelemode
2710 if (scelemode.eq.0) then
2712 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
2714 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
2715 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
2716 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
2717 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
2719 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
2722 ! For the GB potential convert sigma'**2 into chi'
2725 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2729 write (iout,'(/a/)') 'Parameters of the BP potential:'
2730 write (iout,'(a/)') 'The epsilon array:'
2731 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2732 write (iout,'(/a)') 'One-body parameters:'
2733 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
2735 write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),&
2736 chip(i),alp(i),i=1,ntyp)
2739 ! print *,ntyp,"NTYP"
2740 allocate(icharge(ntyp1))
2741 ! print *,ntyp,icharge(i)
2743 read (isidep,*) (icharge(i),i=1,ntyp)
2744 print *,ntyp,icharge(i)
2745 ! if(.not.allocated(eps)) allocate(eps(-ntyp
2746 !c write (2,*) "icharge",(icharge(i),i=1,ntyp)
2747 allocate(alphapol(ntyp,ntyp),epshead(ntyp,ntyp),sig0head(ntyp,ntyp))
2748 allocate(sigiso1(ntyp,ntyp),rborn(ntyp,ntyp),sigmap1(ntyp,ntyp))
2749 allocate(sigmap2(ntyp,ntyp),sigiso2(ntyp,ntyp))
2750 allocate(chis(ntyp,ntyp),wquad(ntyp,ntyp),chipp(ntyp,ntyp))
2751 allocate(epsintab(ntyp,ntyp))
2752 allocate(dtail(2,ntyp,ntyp))
2753 allocate(alphasur(4,ntyp,ntyp),alphiso(4,ntyp,ntyp))
2754 allocate(wqdip(2,ntyp,ntyp))
2755 allocate(wstate(4,ntyp,ntyp))
2756 allocate(dhead(2,2,ntyp,ntyp))
2757 allocate(nstate(ntyp,ntyp))
2758 allocate(debaykap(ntyp,ntyp))
2760 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2761 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2765 ! write (*,*) "Im in ALAB", i, " ", j
2767 eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & !6 w tej linii
2768 (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), & !6 w tej linii
2769 chis(i,j),chis(j,i), & !2 w tej linii
2770 nstate(i,j),(wstate(k,i,j),k=1,4), & !5 w tej lini - 1 integer pierwszy
2771 dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),& ! 4 w tej linii
2772 dtail(1,i,j),dtail(2,i,j), & ! 2 w tej lini
2773 epshead(i,j),sig0head(i,j), & ! 2 w tej linii
2774 rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), & ! 5 w tej linii
2775 alphapol(i,j),alphapol(j,i), & ! 2 w tej linii
2776 (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j) !8 w tej linii
2777 IF ((LaTeX).and.(i.gt.24)) then
2778 write (2,'(2a4,1h&,14(f8.2,1h&),23(f8.2,1h&))') restyp(i,1),restyp(j,1), &
2779 eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & !6 w tej linii
2780 (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), & !6 w tej linii
2781 chis(i,j),chis(j,i) !2 w tej linii
2783 print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i), wqdip(1,i,j)
2787 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
2791 IF ((LaTeX).and.(i.gt.24)) then
2792 write (2,'(2a4,1h&,14(f8.2,1h&),23(f8.2,1h&))') restyp(i,1),restyp(j,1), &
2793 dhead(1,1,i,j),dhead(2,1,i,j),& ! 2 w tej linii
2794 dtail(1,i,j),dtail(2,i,j), & ! 2 w tej lini
2795 epshead(i,j),sig0head(i,j), & ! 2 w tej linii
2796 rborn(i,j),rborn(j,i), & ! 3 w tej linii
2797 alphapol(i,j),alphapol(j,i), & ! 2 w tej linii
2798 (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j) !8 w tej linii
2805 sigma(i,j) = sigma(j,i)
2806 sigmap1(i,j)=sigmap1(j,i)
2807 sigmap2(i,j)=sigmap2(j,i)
2808 sigiso1(i,j)=sigiso1(j,i)
2809 sigiso2(i,j)=sigiso2(j,i)
2810 ! print *,"ATU",sigma(j,i),sigma(i,j),i,j
2811 nstate(i,j) = nstate(j,i)
2812 dtail(1,i,j) = dtail(2,j,i)
2813 dtail(2,i,j) = dtail(1,j,i)
2815 alphasur(k,i,j) = alphasur(k,j,i)
2816 wstate(k,i,j) = wstate(k,j,i)
2817 alphiso(k,i,j) = alphiso(k,j,i)
2820 dhead(2,1,i,j) = dhead(1,1,j,i)
2821 dhead(2,2,i,j) = dhead(1,2,j,i)
2822 dhead(1,1,i,j) = dhead(2,1,j,i)
2823 dhead(1,2,i,j) = dhead(2,2,j,i)
2825 epshead(i,j) = epshead(j,i)
2826 sig0head(i,j) = sig0head(j,i)
2829 wqdip(k,i,j) = wqdip(k,j,i)
2832 wquad(i,j) = wquad(j,i)
2833 epsintab(i,j) = epsintab(j,i)
2834 debaykap(i,j)=debaykap(j,i)
2835 ! if (epsintab(i,j).ne.1.0) print *,"WHAT?",i,j,epsintab(i,j)
2842 !--------------------- GBV potential -----------------------------------
2844 ! 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2845 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2846 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2847 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2849 write (iout,'(/a/)') 'Parameters of the GBV potential:'
2850 write (iout,'(a/)') 'The epsilon array:'
2851 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2852 write (iout,'(/a)') 'One-body parameters:'
2853 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
2854 's||/s_|_^2',' chip ',' alph '
2855 write (iout,'(a3,6x,5f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2856 sigii(i),chip(i),alp(i),i=1,ntyp)
2859 write(iout,*)"Wrong ipot"
2865 !-----------------------------------------------------------------------
2866 ! Calculate the "working" parameters of SC interactions.
2868 !el from module energy - COMMON.INTERACT-------
2869 allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1))
2870 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2871 allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2872 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2873 allocate(r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2874 allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2876 allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2882 if (scelemode.eq.0) then
2891 sc_aa_tube_par(:)=0.0d0
2892 sc_bb_tube_par(:)=0.0d0
2894 !--------------------------------
2899 epslip(i,j)=epslip(j,i)
2902 if (scelemode.eq.0) then
2905 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2906 sigma(j,i)=sigma(i,j)
2907 rs0(i,j)=dwa16*sigma(i,j)
2912 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2913 'Working parameters of the SC interactions:',&
2914 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2919 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2921 ! print *,"SIGMA ZLA?",sigma(i,j)
2929 sigeps=dsign(1.0D0,epsij)
2931 aa_aq(i,j)=epsij*rrij*rrij
2932 ! print *,"ADASKO",epsij,rrij,expon
2933 bb_aq(i,j)=-sigeps*epsij*rrij
2934 aa_aq(j,i)=aa_aq(i,j)
2935 bb_aq(j,i)=bb_aq(i,j)
2936 epsijlip=epslip(i,j)
2937 sigeps=dsign(1.0D0,epsijlip)
2938 epsijlip=dabs(epsijlip)
2939 aa_lip(i,j)=epsijlip*rrij*rrij
2940 bb_lip(i,j)=-sigeps*epsijlip*rrij
2941 aa_lip(j,i)=aa_lip(i,j)
2942 bb_lip(j,i)=bb_lip(i,j)
2943 !C write(iout,*) aa_lip
2944 if ((ipot.gt.2).and. (scelemode.eq.0)) then
2945 sigt1sq=sigma0(i)**2
2946 sigt2sq=sigma0(j)**2
2949 ratsig1=sigt2sq/sigt1sq
2950 ratsig2=1.0D0/ratsig1
2951 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2952 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2953 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2957 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2958 sigmaii(i,j)=rsum_max
2959 sigmaii(j,i)=rsum_max
2961 ! sigmaii(i,j)=r0(i,j)
2962 ! sigmaii(j,i)=r0(i,j)
2964 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2965 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2966 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2967 augm(i,j)=epsij*r_augm**(2*expon)
2968 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2975 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2976 restyp(i,1),restyp(j,1),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2977 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2982 allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2)))
2983 allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2984 allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2985 allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2986 allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2987 allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2988 allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2989 allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2990 allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2991 allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2992 allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2993 allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2994 allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2995 allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2996 allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2997 allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2)
3006 read (isidep_nucl,*) ipot_nucl
3007 ! print *,"TU?!",ipot_nucl
3008 if (ipot_nucl.eq.1) then
3009 do i=1,ntyp_molec(2)
3010 do j=i,ntyp_molec(2)
3011 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),&
3012 elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j)
3016 do i=1,ntyp_molec(2)
3017 do j=i,ntyp_molec(2)
3018 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),&
3019 chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),&
3020 elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j)
3024 ! rpp(1,1)=2**(1.0/6.0)*5.16158
3025 do i=1,ntyp_molec(2)
3026 do j=i,ntyp_molec(2)
3027 rrij=sigma_nucl(i,j)
3031 epsij=4*eps_nucl(i,j)
3032 sigeps=dsign(1.0D0,epsij)
3034 aa_nucl(i,j)=epsij*rrij*rrij
3035 bb_nucl(i,j)=-sigeps*epsij*rrij
3036 ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij)
3037 ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij
3038 ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij
3039 ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij
3040 sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- &
3041 2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i)))
3044 aa_nucl(i,j)=aa_nucl(j,i)
3045 bb_nucl(i,j)=bb_nucl(j,i)
3046 ael3_nucl(i,j)=ael3_nucl(j,i)
3047 ael6_nucl(i,j)=ael6_nucl(j,i)
3048 ael63_nucl(i,j)=ael63_nucl(j,i)
3049 ael32_nucl(i,j)=ael32_nucl(j,i)
3050 elpp3_nucl(i,j)=elpp3_nucl(j,i)
3051 elpp6_nucl(i,j)=elpp6_nucl(j,i)
3052 elpp63_nucl(i,j)=elpp63_nucl(j,i)
3053 elpp32_nucl(i,j)=elpp32_nucl(j,i)
3054 eps_nucl(i,j)=eps_nucl(j,i)
3055 sigma_nucl(i,j)=sigma_nucl(j,i)
3056 sigmaii_nucl(i,j)=sigmaii_nucl(j,i)
3060 write(iout,*) "tube param"
3061 read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
3062 ccavtubpep,dcavtubpep,tubetranenepep
3063 sigmapeptube=sigmapeptube**6
3064 sigeps=dsign(1.0D0,epspeptube)
3065 epspeptube=dabs(epspeptube)
3066 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
3067 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
3068 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
3070 read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
3071 ccavtub(i),dcavtub(i),tubetranene(i)
3072 sigmasctube=sigmasctube**6
3073 sigeps=dsign(1.0D0,epssctube)
3074 epssctube=dabs(epssctube)
3075 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
3076 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
3077 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
3079 !-----------------READING SC BASE POTENTIALS-----------------------------
3080 allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))
3081 allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
3082 allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
3083 allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
3084 allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
3085 allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
3086 allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
3087 allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
3088 allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
3089 allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
3090 allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
3091 allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
3092 allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
3093 allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
3094 allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
3095 allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
3097 write (iout,*) "ESCBASEPARM"
3098 do i=1,ntyp_molec(1)
3099 do j=1,ntyp_molec(2) ! without U then we will take T for U
3100 ! write (*,*) "Im in ", i, " ", j
3101 read(isidep_scbase,*) &
3102 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
3103 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
3104 ! write(*,*) "eps",eps_scbase(i,j)
3105 read(isidep_scbase,*) &
3106 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
3107 chis_scbase(i,j,1),chis_scbase(i,j,2)
3108 read(isidep_scbase,*) &
3109 dhead_scbasei(i,j), &
3110 dhead_scbasej(i,j), &
3111 rborn_scbasei(i,j),rborn_scbasej(i,j)
3112 read(isidep_scbase,*) &
3113 (wdipdip_scbase(k,i,j),k=1,3), &
3114 (wqdip_scbase(k,i,j),k=1,2)
3115 read(isidep_scbase,*) &
3116 alphapol_scbase(i,j), &
3117 epsintab_scbase(i,j)
3118 if (chi_scbase(i,j,2).gt.0.9) chi_scbase(i,j,2)=0.9
3119 if (chi_scbase(i,j,1).gt.0.9) chi_scbase(i,j,1)=0.9
3120 if (chipp_scbase(i,j,2).gt.0.9) chipp_scbase(i,j,2)=0.9
3121 if (chipp_scbase(i,j,1).gt.0.9) chipp_scbase(i,j,1)=0.9
3122 if (chi_scbase(i,j,2).lt.-0.9) chi_scbase(i,j,2)=-0.9
3123 if (chi_scbase(i,j,1).lt.-0.9) chi_scbase(i,j,1)=-0.9
3124 if (chipp_scbase(i,j,2).lt.-0.9) chipp_scbase(i,j,2)=-0.9
3125 if (chipp_scbase(i,j,1).lt.-0.9) chipp_scbase(i,j,1)=-0.9
3127 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
3128 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
3129 write(*,*) "eps",eps_scbase(i,j)
3131 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
3132 chis_scbase(i,j,1),chis_scbase(i,j,2)
3134 dhead_scbasei(i,j), &
3135 dhead_scbasej(i,j), &
3136 rborn_scbasei(i,j),rborn_scbasej(i,j)
3138 (wdipdip_scbase(k,i,j),k=1,3), &
3139 (wqdip_scbase(k,i,j),k=1,2)
3141 alphapol_scbase(i,j), &
3142 epsintab_scbase(i,j)
3147 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
3148 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
3149 write(*,*) "eps",eps_scbase(i,j)
3151 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
3152 chis_scbase(i,j,1),chis_scbase(i,j,2)
3154 dhead_scbasei(i,j), &
3155 dhead_scbasej(i,j), &
3156 rborn_scbasei(i,j),rborn_scbasej(i,j)
3158 (wdipdip_scbase(k,i,j),k=1,3), &
3159 (wqdip_scbase(k,i,j),k=1,2)
3161 alphapol_scbase(i,j), &
3162 epsintab_scbase(i,j)
3165 allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
3166 allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
3168 do i=1,ntyp_molec(1)
3169 do j=1,ntyp_molec(2)-1
3170 epsij=eps_scbase(i,j)
3171 rrij=sigma_scbase(i,j)
3176 sigeps=dsign(1.0D0,epsij)
3178 aa_scbase(i,j)=epsij*rrij*rrij
3179 bb_scbase(i,j)=-sigeps*epsij*rrij
3182 !-----------------READING PEP BASE POTENTIALS-------------------
3183 allocate(eps_pepbase(ntyp_molec(2)))
3184 allocate(sigma_pepbase(ntyp_molec(2)))
3185 allocate(chi_pepbase(ntyp_molec(2),2))
3186 allocate(chipp_pepbase(ntyp_molec(2),2))
3187 allocate(alphasur_pepbase(4,ntyp_molec(2)))
3188 allocate(sigmap1_pepbase(ntyp_molec(2)))
3189 allocate(sigmap2_pepbase(ntyp_molec(2)))
3190 allocate(chis_pepbase(ntyp_molec(2),2))
3191 allocate(wdipdip_pepbase(3,ntyp_molec(2)))
3193 write (iout,*) "EPEPBASEPARM"
3195 do j=1,ntyp_molec(2) ! without U then we will take T for U
3196 write (*,*) "Im in ", i, " ", j
3197 read(isidep_pepbase,*) &
3198 eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
3199 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
3200 if (chi_pepbase(j,2).gt.0.9) chi_pepbase(j,2)=0.9
3201 if (chi_pepbase(j,1).gt.0.9) chi_pepbase(j,1)=0.9
3202 if (chipp_pepbase(j,2).gt.0.9) chipp_pepbase(j,2)=0.9
3203 if (chipp_pepbase(j,1).gt.0.9) chipp_pepbase(j,1)=0.9
3204 if (chi_pepbase(j,2).lt.-0.9) chi_pepbase(j,2)=-0.9
3205 if (chi_pepbase(j,1).lt.-0.9) chi_pepbase(j,1)=-0.9
3206 if (chipp_pepbase(j,2).lt.-0.9) chipp_pepbase(j,2)=-0.9
3207 if (chipp_pepbase(j,1).lt.-0.9) chipp_pepbase(j,1)=-0.9
3209 write(*,*) "eps",eps_pepbase(j)
3210 read(isidep_pepbase,*) &
3211 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
3212 chis_pepbase(j,1),chis_pepbase(j,2)
3213 read(isidep_pepbase,*) &
3214 (wdipdip_pepbase(k,j),k=1,3)
3215 write(iout,*) eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
3216 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
3218 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
3219 chis_pepbase(j,1),chis_pepbase(j,2)
3221 (wdipdip_pepbase(k,j),k=1,3)
3225 write(iout,*) eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
3226 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
3228 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
3229 chis_pepbase(j,1),chis_pepbase(j,2)
3231 (wdipdip_pepbase(k,j),k=1,3)
3233 allocate(aa_pepbase(ntyp_molec(2)))
3234 allocate(bb_pepbase(ntyp_molec(2)))
3236 do j=1,ntyp_molec(2)-1
3237 epsij=eps_pepbase(j)
3238 rrij=sigma_pepbase(j)
3243 sigeps=dsign(1.0D0,epsij)
3245 aa_pepbase(j)=epsij*rrij*rrij
3246 bb_pepbase(j)=-sigeps*epsij*rrij
3248 !--------------READING SC PHOSPHATE-------------------------------------
3249 allocate(eps_scpho(ntyp_molec(1)))
3250 allocate(sigma_scpho(ntyp_molec(1)))
3251 allocate(chi_scpho(ntyp_molec(1),2))
3252 allocate(chipp_scpho(ntyp_molec(1),2))
3253 allocate(alphasur_scpho(4,ntyp_molec(1)))
3254 allocate(sigmap1_scpho(ntyp_molec(1)))
3255 allocate(sigmap2_scpho(ntyp_molec(1)))
3256 allocate(chis_scpho(ntyp_molec(1),2))
3257 allocate(wqq_scpho(ntyp_molec(1)))
3258 allocate(wqdip_scpho(2,ntyp_molec(1)))
3259 allocate(alphapol_scpho(ntyp_molec(1)))
3260 allocate(epsintab_scpho(ntyp_molec(1)))
3261 allocate(dhead_scphoi(ntyp_molec(1)))
3262 allocate(rborn_scphoi(ntyp_molec(1)))
3263 allocate(rborn_scphoj(ntyp_molec(1)))
3264 allocate(alphi_scpho(ntyp_molec(1)))
3268 do j=1,ntyp_molec(1) ! without U then we will take T for U
3269 write (*,*) "Im in scpho ", i, " ", j
3270 read(isidep_scpho,*) &
3271 eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
3272 chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
3273 write(*,*) "eps",eps_scpho(j)
3274 read(isidep_scpho,*) &
3275 (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
3276 chis_scpho(j,1),chis_scpho(j,2)
3277 read(isidep_scpho,*) &
3278 (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
3279 read(isidep_scpho,*) &
3280 epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), &
3282 if (chi_scpho(j,2).gt.0.9) chi_scpho(j,2)=0.9
3283 if (chi_scpho(j,1).gt.0.9) chi_scpho(j,1)=0.9
3284 if (chipp_scpho(j,2).gt.0.9) chipp_scpho(j,2)=0.9
3285 if (chipp_scpho(j,1).gt.0.9) chipp_scpho(j,1)=0.9
3286 if (chi_scpho(j,2).lt.-0.9) chi_scpho(j,2)=-0.9
3287 if (chi_scpho(j,1).lt.-0.9) chi_scpho(j,1)=-0.9
3288 if (chipp_scpho(j,2).lt.-0.9) chipp_scpho(j,2)=-0.9
3289 if (chipp_scpho(j,1).lt.-0.9) chipp_scpho(j,1)=-0.9
3293 allocate(aa_scpho(ntyp_molec(1)))
3294 allocate(bb_scpho(ntyp_molec(1)))
3296 do j=1,ntyp_molec(1)
3303 sigeps=dsign(1.0D0,epsij)
3305 aa_scpho(j)=epsij*rrij*rrij
3306 bb_scpho(j)=-sigeps*epsij*rrij
3310 read(isidep_peppho,*) &
3311 eps_peppho,sigma_peppho
3312 read(isidep_peppho,*) &
3313 (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
3314 read(isidep_peppho,*) &
3315 (wqdip_peppho(k),k=1,2)
3323 sigeps=dsign(1.0D0,epsij)
3325 aa_peppho=epsij*rrij*rrij
3326 bb_peppho=-sigeps*epsij*rrij
3329 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
3334 ! Define the SC-p interaction constants (hard-coded; old style)
3337 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
3339 ! aad(i,1)=0.3D0*4.0D0**12
3340 ! Following line for constants currently implemented
3341 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
3342 aad(i,1)=1.5D0*4.0D0**12
3343 ! aad(i,1)=0.17D0*5.6D0**12
3345 ! "Soft" SC-p repulsion
3347 ! Following line for constants currently implemented
3348 ! aad(i,1)=0.3D0*4.0D0**6
3349 ! "Hard" SC-p repulsion
3350 bad(i,1)=3.0D0*4.0D0**6
3351 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
3360 ! 8/9/01 Read the SC-p interaction constants from file
3363 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
3366 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
3367 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
3368 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
3369 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
3373 write (iout,*) "Parameters of SC-p interactions:"
3375 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
3376 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
3381 allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
3383 do i=1,ntyp_molec(2)
3384 read (iscpp_nucl,*,end=118,err=118) eps_scp_nucl(i),rscp_nucl(i)
3386 do i=1,ntyp_molec(2)
3387 aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
3388 bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
3390 r0pp=1.12246204830937298142*5.16158
3396 ! Define the constants of the disulfide bridge
3400 ! Old arbitrary potential - commented out.
3405 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
3406 ! energy surface of diethyl disulfide.
3407 ! A. Liwo and U. Kozlowska, 11/24/03
3423 allocate(ichargelipid(ntyp_molec(4)))
3424 allocate(lip_angle_force(ntyp_molec(4),ntyp_molec(4),ntyp_molec(4)))
3425 allocate(lip_angle_angle(ntyp_molec(4),ntyp_molec(4),ntyp_molec(4)))
3426 allocate(lip_bond(ntyp_molec(4),ntyp_molec(4)))
3427 allocate(lip_eps(ntyp_molec(4),ntyp_molec(4)))
3428 allocate(lip_sig(ntyp_molec(4),ntyp_molec(4)))
3429 kjtokcal=0.2390057361
3431 !HERE THE MASS of MARTINI
3432 write(*,*) "before MARTINI PARAM"
3433 do i=1,ntyp_molec(4)
3439 !relative dielectric constant = 15 for implicit screening
3440 k_coulomb_lip=332.0d0/15.0d0
3441 !kbond = 1250 kJ/(mol*nm*2)
3442 kbondlip=1250.0d0*kjtokcal/100.0d0
3444 lip_angle_force=0.0d0
3445 if (DRY_MARTINI.gt.0) then
3446 lip_angle_force(3,12,12)=35.0*kjtokcal!*krad2
3447 lip_angle_force(3,12,18)=35.0*kjtokcal!*krad2
3448 lip_angle_force(3,18,16)=35.0*kjtokcal!*krad2
3449 lip_angle_force(12,18,16)=35.0*kjtokcal!*krad2
3450 lip_angle_force(18,16,18)=45.0*kjtokcal!*krad2
3451 lip_angle_force(16,18,18)=35.0*kjtokcal!*krad2
3452 lip_angle_force(12,18,18)=35.0*kjtokcal!*krad2
3453 lip_angle_force(18,18,18)=35.0*kjtokcal!*krad2
3455 lip_angle_force(3,12,12)=25.0*kjtokcal!*krad2
3456 lip_angle_force(3,12,18)=25.0*kjtokcal!*krad2
3457 lip_angle_force(3,18,16)=25.0*kjtokcal!*krad2
3458 lip_angle_force(12,18,16)=25.0*kjtokcal!*krad2
3459 lip_angle_force(18,16,18)=45.0*kjtokcal!*krad2
3460 lip_angle_force(16,18,18)=25.0*kjtokcal!*krad2
3461 lip_angle_force(12,18,18)=25.0*kjtokcal!*krad2
3462 lip_angle_force(18,18,18)=25.0*kjtokcal!*krad2
3464 lip_angle_angle=0.0d0
3465 lip_angle_angle(3,12,12)=120.0/krad
3466 lip_angle_angle(3,12,18)=180.0/krad
3467 lip_angle_angle(3,18,16)=180.0/krad
3468 lip_angle_angle(12,18,16)=180.0/krad
3469 lip_angle_angle(18,16,18)=120.0/krad
3470 lip_angle_angle(16,18,18)=180.0/krad
3471 lip_angle_angle(12,18,18)=180.0/krad
3472 lip_angle_angle(18,18,18)=180.0/krad
3473 read(ilipbond,*) temp1
3475 read(ilipbond,*) temp1, lip_bond(i,1), &
3476 lip_bond(i,2),lip_bond(i,3),lip_bond(i,4),lip_bond(i,5), &
3477 lip_bond(i,6),lip_bond(i,7),lip_bond(i,8),lip_bond(i,9), &
3478 lip_bond(i,10),lip_bond(i,11),lip_bond(i,12),lip_bond(i,13), &
3479 lip_bond(i,14),lip_bond(i,15),lip_bond(i,16),lip_bond(i,17), &
3482 lip_bond(i,j)=lip_bond(i,j)*10
3486 read(ilipnonbond,*) (ichargelipid(i),i=1,ntyp_molec(4))
3487 read(ilipnonbond,*) temp1
3489 read(ilipnonbond,*) temp1, lip_eps(i,1), &
3490 lip_eps(i,2),lip_eps(i,3),lip_eps(i,4),lip_eps(i,5), &
3491 lip_eps(i,6),lip_eps(i,7),lip_eps(i,8),lip_eps(i,9), &
3492 lip_eps(i,10),lip_eps(i,11),lip_eps(i,12),lip_eps(i,13), &
3493 lip_eps(i,14),lip_eps(i,15),lip_eps(i,16),lip_eps(i,17), &
3495 ! write(*,*) i, lip_eps(i,18)
3497 lip_eps(i,j)=lip_eps(i,j)*kjtokcal
3500 read(ilipnonbond,*) temp1
3502 read(ilipnonbond,*) temp1,lip_sig(i,1), &
3503 lip_sig(i,2),lip_sig(i,3),lip_sig(i,4),lip_sig(i,5), &
3504 lip_sig(i,6),lip_sig(i,7),lip_sig(i,8),lip_sig(i,9), &
3505 lip_sig(i,10),lip_sig(i,11),lip_sig(i,12),lip_sig(i,13), &
3506 lip_sig(i,14),lip_sig(i,15),lip_sig(i,16),lip_sig(i,17), &
3509 lip_sig(i,j)=lip_sig(i,j)*10.0
3512 write(*,*) "after MARTINI PARAM"
3516 allocate(alphapolcat(ntyp,ntyp),epsheadcat(ntyp,ntyp),sig0headcat(ntyp,ntyp))
3517 allocate(alphapolcat2(ntyp,ntyp))
3518 allocate(sigiso1cat(ntyp,ntyp),rborn1cat(ntyp,ntyp),rborn2cat(ntyp,ntyp),sigmap1cat(ntyp,ntyp))
3519 allocate(sigmap2cat(ntyp,ntyp),sigiso2cat(ntyp,ntyp))
3520 allocate(chis1cat(ntyp,ntyp),chis2cat(ntyp,ntyp),wquadcat(ntyp,ntyp),chipp1cat(ntyp,ntyp),chipp2cat(ntyp,ntyp))
3521 allocate(epsintabcat(ntyp,ntyp))
3522 allocate(dtailcat(2,ntyp,ntyp))
3523 allocate(alphasurcat(4,ntyp,ntyp),alphisocat(4,ntyp,ntyp))
3524 allocate(wqdipcat(2,ntyp,ntyp))
3525 allocate(wstatecat(4,ntyp,ntyp))
3526 allocate(dheadcat(2,2,ntyp,ntyp))
3527 allocate(nstatecat(ntyp,ntyp))
3528 allocate(debaykapcat(ntyp,ntyp))
3530 if (.not.allocated(epscat)) allocate (epscat(0:ntyp1,0:ntyp1))
3531 if (.not.allocated(sigmacat)) allocate(sigmacat(0:ntyp1,0:ntyp1))
3532 ! if (.not.allocated(chicat)) allocate(chicat(ntyp1,ntyp1)) !(ntyp,ntyp)
3533 if (.not.allocated(chi1cat)) allocate(chi1cat(ntyp1,ntyp1)) !(ntyp,ntyp)
3534 if (.not.allocated(chi2cat)) allocate(chi2cat(ntyp1,ntyp1)) !(ntyp,ntyp)
3537 if (.not.allocated(ichargecat))&
3538 allocate (ichargecat(-ntyp_molec(5):ntyp_molec(5)))
3539 write(*,*) "before ions",oldion
3542 ! i to SC, j to jon, isideocat - nazwa pliku z ktorego czytam parametry
3543 if (oldion.eq.0) then
3544 if (.not.allocated(icharge)) then ! this mean you are oprating in old sc-sc mode
3545 allocate(icharge(1:ntyp1))
3546 read(iion,*) (icharge(i),i=1,ntyp)
3550 print *,ntyp_molec(5)
3551 do i=-ntyp_molec(5),ntyp_molec(5)
3552 read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
3553 print *,msc(i,5),restok(i,5)
3559 do j=1,ntyp_molec(5)-1 ! this is without Zn will be modified for ALL tranistion metals
3561 ! do j=1,ntyp_molec(5)
3562 ! write (*,*) "Im in ALAB", i, " ", j
3564 epscat(i,j),sigmacat(i,j), &
3565 ! chicat(i,j),chicat(j,i),chippcat(i,j),chippcat(j,i), &
3566 chi1cat(i,j),chi2cat(i,j),chipp1cat(i,j),chipp2cat(i,j), &
3568 (alphasurcat(k,i,j),k=1,4),sigmap1cat(i,j),sigmap2cat(i,j),&
3569 ! chiscat(i,j),chiscat(j,i), &
3570 chis1cat(i,j),chis2cat(i,j), &
3572 nstatecat(i,j),(wstatecat(k,i,j),k=1,4), & !5 w tej lini - 1 integer pierwszy
3573 dheadcat(1,1,i,j),dheadcat(1,2,i,j),dheadcat(2,1,i,j),dheadcat(2,2,i,j),&
3574 dtailcat(1,i,j),dtailcat(2,i,j), &
3575 epsheadcat(i,j),sig0headcat(i,j), &
3577 ! rborncat(i,j),rborncat(j,i),&
3578 rborn1cat(i,j),rborn2cat(i,j),&
3579 (wqdipcat(k,i,j),k=1,2), &
3580 alphapolcat(i,j),alphapolcat2(j,i), &
3581 (alphisocat(k,i,j),k=1,4),sigiso1cat(i,j),sigiso2cat(i,j),epsintabcat(i,j),debaykapcat(i,j)
3583 if (chi1cat(i,j).gt.0.9) write (*,*) "WTF ANISO", i,j, chi1cat(i,j)
3584 ! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
3586 ! write (iout,*) 'i= ', i, ' j= ', j
3587 ! write (iout,*) 'epsi0= ', epscat(1,j)
3588 ! write (iout,*) 'sigma0= ', sigmacat(1,j)
3589 ! write (iout,*) 'chi(i,j)= ', chicat(1,j)
3590 ! write (iout,*) 'chi(j,i)= ', chicat(j,1)
3591 ! write (iout,*) 'chip(i,j)= ', chippcat(1,j)
3592 ! write (iout,*) 'chip(j,i)= ', chippcat(j,1)
3593 ! write (iout,*) 'alphasur1= ', alphasurcat(1,1,j)
3594 ! write (iout,*) 'alphasur2= ', alphasurcat(2,1,j)
3595 ! write (iout,*) 'alphasur3= ', alphasurcat(3,1,j)
3596 ! write (iout,*) 'alphasur4= ', alphasurcat(4,1,j)
3597 ! write (iout,*) 'sig1= ', sigmap1cat(1,j)
3598 ! write (iout,*) 'chis(i,j)= ', chiscat(1,j)
3599 ! write (iout,*) 'chis(j,i)= ', chiscat(j,1)
3600 ! write (iout,*) 'dhead= ', dheadcat(1,1,1,j)
3601 ! write (iout,*) 'a1= ', rborncat(j,1)
3602 ! write (iout,*) 'a2= ', rborncat(1,j)
3603 ! write (iout,*) 'epsin= ', epsintabcat(1,j), epsintabcat(j,1)
3604 ! write (iout,*) 'alphapol1= ', alphapolcat(1,j)
3605 ! write (iout,*) 'w1= ', wqdipcat(1,1,j)
3606 ! write (iout,*) 'w2= ', wqdipcat(2,1,j)
3610 ! If ((i.eq.1).and.(j.eq.27)) then
3611 ! write (iout,*) 'SEP'
3612 ! Write (iout,*) 'w1= ', wqdipcat(1,1,27)
3613 ! Write (iout,*) 'w2= ', wqdipcat(2,1,27)
3618 allocate(aa_aq_cat(-ntyp:ntyp,ntyp),bb_aq_cat(-ntyp:ntyp,ntyp))
3620 do j=1,ntyp_molec(5)
3624 sigeps=dsign(1.0D0,epsij)
3626 aa_aq_cat(i,j)=epsij*rrij*rrij
3627 bb_aq_cat(i,j)=-sigeps*epsij*rrij
3632 do j=1,ntyp_molec(5)-1
3634 write (iout,*) 'i= ', i, ' j= ', j
3635 write (iout,*) 'epsi0= ', epscat(i,j)
3636 write (iout,*) 'sigma0= ', sigmacat(i,j)
3637 write (iout,*) 'chi1= ', chi1cat(i,j)
3638 write (iout,*) 'chi1= ', chi2cat(i,j)
3639 write (iout,*) 'chip1= ', chipp1cat(i,j)
3640 write (iout,*) 'chip2= ', chipp2cat(i,j)
3641 write (iout,*) 'alphasur1= ', alphasurcat(1,i,j)
3642 write (iout,*) 'alphasur2= ', alphasurcat(2,i,j)
3643 write (iout,*) 'alphasur3= ', alphasurcat(3,i,j)
3644 write (iout,*) 'alphasur4= ', alphasurcat(4,i,j)
3645 write (iout,*) 'sig1= ', sigmap1cat(i,j)
3646 write (iout,*) 'sig2= ', sigmap2cat(i,j)
3647 write (iout,*) 'chis1= ', chis1cat(i,j)
3648 write (iout,*) 'chis1= ', chis2cat(i,j)
3649 write (iout,*) 'nstatecat(i,j)= ', nstatecat(i,j)
3650 write (iout,*) 'wstatecat(k,i,j),k=1= ',wstatecat(1,i,j)
3651 write (iout,*) 'dhead= ', dheadcat(1,1,i,j)
3652 write (iout,*) 'dhead2= ', dheadcat(1,2,i,j)
3653 write (iout,*) 'a1= ', rborn1cat(i,j)
3654 write (iout,*) 'a2= ', rborn2cat(i,j)
3655 write (iout,*) 'epsin= ', epsintabcat(i,j), epsintabcat(j,i)
3656 write (iout,*) 'alphapol1= ', alphapolcat(i,j)
3657 write (iout,*) 'alphapol2= ', alphapolcat2(i,j)
3658 write (iout,*) 'w1= ', wqdipcat(1,i,j)
3659 write (iout,*) 'w2= ', wqdipcat(2,i,j)
3660 write (iout,*) 'debaykapcat(i,j)= ', debaykapcat(i,j)
3663 If ((i.eq.1).and.(j.eq.27)) then
3664 write (iout,*) 'SEP'
3665 Write (iout,*) 'w1= ', wqdipcat(1,1,27)
3666 Write (iout,*) 'w2= ', wqdipcat(2,1,27)
3673 ! read number of Zn2+
3674 ! here two denotes the Zn2+ and Cu2+
3675 write(iout,*) "before TRANPARM"
3676 allocate(aomicattr(0:3,2))
3677 allocate(athetacattran(0:6,5,2))
3678 allocate(agamacattran(3,5,2))
3679 allocate(acatshiftdsc(5,2))
3680 allocate(bcatshiftdsc(5,2))
3681 allocate(demorsecat(5,2))
3682 allocate(alphamorsecat(5,2))
3683 allocate(x0catleft(5,2))
3684 allocate(x0catright(5,2))
3685 allocate(x0cattrans(5,2))
3686 allocate(ntrantyp(2))
3687 do i=1,1 ! currently only Zn2+
3689 read(iiontran,*) ntrantyp(i)
3692 !ASP| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -.5
3693 !CYS| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0left x0right x0transi
3694 !GLU| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -0.5
3695 !HIS| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -.5
3696 read(iiontran,*) (aomicattr(j,i),j=0,3)
3698 read (iiontran,*,err=123,end=123) (agamacattran(k,j,i),k=1,3),&
3699 (athetacattran(k,j,i),k=0,6),acatshiftdsc(j,i),bcatshiftdsc(j,i),&
3700 demorsecat(j,i),alphamorsecat(j,i),x0catleft(j,i),x0catright(j,i),&
3705 write (iout,'(/a)') "Disulfide bridge parameters:"
3706 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
3707 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
3708 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
3709 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
3712 if (shield_mode.gt.0) then
3713 pi=4.0D0*datan(1.0D0)
3714 !C VSolvSphere the volume of solving sphere
3716 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3717 !C there will be no distinction between proline peptide group and normal peptide
3718 !C group in case of shielding parameters
3719 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
3720 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
3721 write (iout,*) VSolvSphere,VSolvSphere_div
3722 !C long axis of side chain
3724 long_r_sidechain(i)=vbldsc0(1,i)
3725 ! if (scelemode.eq.0) then
3726 short_r_sidechain(i)=sigma(i,i)/sqrt(2.0)
3727 if (short_r_sidechain(i).eq.0.0) short_r_sidechain(i)=0.2
3729 ! short_r_sidechain(i)=sigma(i,i)
3731 write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3738 111 write (iout,*) "Error reading bending energy parameters."
3740 112 write (iout,*) "Error reading rotamer energy parameters."
3742 113 write (iout,*) "Error reading torsional energy parameters."
3744 114 write (iout,*) "Error reading double torsional energy parameters."
3746 115 write (iout,*) &
3747 "Error reading cumulant (multibody energy) parameters."
3749 116 write (iout,*) "Error reading electrostatic energy parameters."
3751 117 write (iout,*) "Error reading side chain interaction parameters."
3753 118 write (iout,*) "Error reading SCp interaction parameters."
3755 119 write (iout,*) "Error reading SCCOR parameters"
3757 121 write (iout,*) "Error in Czybyshev parameters"
3759 123 write(iout,*) "Error in transition metal parameters"
3762 call MPI_Finalize(Ierror)
3766 end subroutine parmread
3768 !-----------------------------------------------------------------------------
3770 !-----------------------------------------------------------------------------
3771 subroutine printmat(ldim,m,n,iout,key,a)
3774 character(len=3),dimension(n) :: key
3775 real(kind=8),dimension(ldim,n) :: a
3777 integer :: i,j,k,m,iout,nlim
3781 write (iout,1000) (key(k),k=i,nlim)
3783 1000 format (/5x,8(6x,a3))
3784 1020 format (/80(1h-)/)
3786 write (iout,1010) key(j),(a(j,k),k=i,nlim)
3789 1010 format (a3,2x,8(f9.4))
3791 end subroutine printmat
3792 !-----------------------------------------------------------------------------
3794 !-----------------------------------------------------------------------------
3796 ! Read the PDB file and convert the peptide geometry into virtual-chain
3799 use energy_data, only: itype,istype
3803 ! use control, only: rescode,sugarcode
3804 ! implicit real*8 (a-h,o-z)
3805 ! include 'DIMENSIONS'
3806 ! include 'COMMON.LOCAL'
3807 ! include 'COMMON.VAR'
3808 ! include 'COMMON.CHAIN'
3809 ! include 'COMMON.INTERACT'
3810 ! include 'COMMON.IOUNITS'
3811 ! include 'COMMON.GEO'
3812 ! include 'COMMON.NAMES'
3813 ! include 'COMMON.CONTROL'
3814 ! include 'COMMON.DISTFIT'
3815 ! include 'COMMON.SETUP'
3816 integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,&
3818 logical :: lprn=.true.,fail
3819 real(kind=8),dimension(3) :: e1,e2,e3
3820 real(kind=8) :: dcj,efree_temp
3821 character(len=3) :: seq,res,res2
3822 character(len=5) :: atom
3823 character(len=80) :: card
3824 real(kind=8),dimension(3,40) :: sccor
3825 integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode,
3826 integer :: isugar,molecprev,firstion
3827 character*1 :: sugar
3829 real(kind=8),dimension(3) :: ccc
3831 integer,dimension(2,maxres/3) :: hfrag_alloc
3832 integer,dimension(4,maxres/3) :: bfrag_alloc
3833 real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
3834 real(kind=8),dimension(:,:), allocatable :: c_temporary
3835 integer,dimension(:,:) , allocatable :: itype_temporary
3836 integer,dimension(:),allocatable :: istype_temp
3843 ! write (2,*) "UNRES_PDB",unres_pdb
3863 !-----------------------------
3864 allocate(hfrag(2,maxres/3)) !(2,maxres/3)
3865 allocate(bfrag(4,maxres/3)) !(4,maxres/3)
3866 if(.not. allocated(istype)) allocate(istype(maxres))
3868 read (ipdbin,'(a80)',end=10) card
3869 write (iout,'(a)') card
3870 if (card(:5).eq.'HELIX') then
3873 read(card(22:25),*) hfrag(1,nhfrag)
3874 read(card(34:37),*) hfrag(2,nhfrag)
3876 if (card(:5).eq.'SHEET') then
3879 read(card(24:26),*) bfrag(1,nbfrag)
3880 read(card(35:37),*) bfrag(2,nbfrag)
3881 !rc----------------------------------------
3882 !rc to be corrected !!!
3883 bfrag(3,nbfrag)=bfrag(1,nbfrag)
3884 bfrag(4,nbfrag)=bfrag(2,nbfrag)
3885 !rc----------------------------------------
3887 if (card(:3).eq.'END') then
3889 else if (card(:3).eq.'TER') then
3894 itype(ires_old,molecule)=ntyp1_molec(molecule)
3895 itype(ires_old-1,molecule)=ntyp1_molec(molecule)
3896 nres_molec(molecule)=nres_molec(molecule)+2
3897 ! if (molecule.eq.4) ires=ires+2
3899 ! write (iout,*) "Chain ended",ires,ishift,ires_old
3902 dc(j,ires)=sccor(j,iii)
3905 call sccenter(ires,iii,sccor)
3909 else if (card(:3).eq.'BRA') then
3913 itype(ires,molecule)=ntyp1_molec(molecule)-1
3914 nres_molec(molecule)=nres_molec(molecule)+1
3918 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
3919 ! Fish out the ATOM cards.
3920 ! write(iout,*) 'card',card(1:20)
3921 ! print *,"ATU ",card(1:6), CARD(3:6)
3923 if (index(card(1:4),'ATOM').gt.0) then
3924 read (card(12:16),*) atom
3925 ! write (iout,*) "! ",atom," !",ires
3926 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
3927 if (card(14:16).eq.'LIP') then
3932 nres_molec(molecule)=nres_molec(molecule)+1
3933 itype(ires,molecule)=ntyp1_molec(molecule)
3942 nres_molec(molecule)=nres_molec(molecule)+1
3943 read (card(18:20),'(a3)') res
3945 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3946 if (UNRES_PDB) then!
3947 itype(ires,molecule)=rescode(ires,res,0,molecule)
3949 itype(ires,molecule)=rescode_lip(res,ires)
3952 read (card(23:26),*) ires
3953 read (card(18:20),'(a3)') res
3954 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
3955 ! & " ires_old",ires_old
3956 ! write (iout,*) "ishift",ishift," ishift1",ishift1
3957 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
3958 if (ires-ishift+ishift1.ne.ires_old) then
3959 ! Calculate the CM of the preceding residue.
3960 ! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
3962 ! write (iout,*) "Calculating sidechain center iii",iii
3965 dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
3968 call sccenter(ires_old,iii,sccor)
3972 ! Start new residue.
3973 if (res.eq.'Cl-' .or. res.eq.'Na+') then
3976 else if (ibeg.eq.1) then
3977 write (iout,*) "BEG ires",ires
3979 if (res.ne.'GLY' .and. res.ne. 'ACE') then
3982 nres_molec(molecule)=nres_molec(molecule)+1
3984 ires=ires-ishift+ishift1
3986 ! write (iout,*) "ishift",ishift," ires",ires,&
3987 ! " ires_old",ires_old
3989 else if (ibeg.eq.2) then
3991 ishift=-ires_old+ires-1 !!!!!
3992 ishift1=ishift1-1 !!!!!
3993 ! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
3994 ires=ires-ishift+ishift1
3995 ! print *,ires,ishift,ishift1
3999 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
4000 ires=ires-ishift+ishift1
4003 ! print *,'atom',ires,atom
4004 if (res.eq.'ACE' .or. res.eq.'NHE') then
4007 if (atom.eq.'CA '.or.atom.eq.'N ') then
4009 itype(ires,molecule)=rescode(ires,res,0,molecule)
4011 ! nres_molec(molecule)=nres_molec(molecule)+1
4015 itype(ires,molecule)=rescode(ires,res2,0,molecule)
4016 ! nres_molec(molecule)=nres_molec(molecule)+1
4017 read (card(19:19),'(a1)') sugar
4018 isugar=sugarcode(sugar,ires)
4019 ! if (ibeg.eq.1) then
4023 ! print *,"ires=",ires,istype(ires)
4029 ires=ires-ishift+ishift1
4031 ! write (iout,*) "ires_old",ires_old," ires",ires
4032 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
4035 ! write (2,*) "ires",ires," res ",res!," ity"!,ity
4036 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
4037 res.eq.'NHE'.and.atom(:2).eq.'HN') then
4038 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
4039 ! print *,ires,ishift,ishift1
4040 ! write (iout,*) "backbone ",atom
4042 write (iout,'(2i3,2x,a,3f8.3)') &
4043 ires,itype(ires,1),res,(c(j,ires),j=1,3)
4046 nres_molec(molecule)=nres_molec(molecule)+1
4048 sccor(j,iii)=c(j,ires)
4050 else if (.not.unres_pdb .and. (atom.eq."C1'" .or. &
4051 atom.eq."C2'" .or. atom.eq."C3'" &
4052 .or. atom.eq."C4'" .or. atom.eq."O4'")) then
4053 read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
4054 !c write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
4055 ! print *,ires,ishift,ishift1
4059 ! sccor(j,iii)=c(j,ires)
4062 c(j,ires)=c(j,ires)+ccc(j)/5.0
4064 print *,counter,molecule
4065 if (counter.eq.5) then
4067 nres_molec(molecule)=nres_molec(molecule)+1
4070 ! sccor(j,iii)=c(j,ires)
4074 ! print *, "ATOM",atom(1:3)
4075 else if (atom.eq."C5'") then
4076 read (card(19:19),'(a1)') sugar
4077 isugar=sugarcode(sugar,ires)
4082 ! print *,ires,istype(ires)
4086 ! print *,"nres_molec(molecule)",nres_molec(molecule),ires
4087 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
4088 nres_molec(molecule)=nres_molec(molecule)+1
4089 print *,"nres_molec(molecule)",nres_molec(molecule),ires
4093 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
4095 else if ((atom.eq."C1'").and.unres_pdb) then
4097 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
4098 ! write (*,*) card(23:27),ires,itype(ires,1)
4099 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
4100 atom.ne.'N' .and. atom.ne.'C' .and. &
4101 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
4102 atom.ne.'OXT' .and. atom(:2).ne.'3H' &
4103 .and. atom.ne.'P '.and. &
4104 atom(1:1).ne.'H' .and. &
4105 atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'&
4107 ! write (iout,*) "sidechain ",atom
4108 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
4109 if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
4110 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
4112 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
4116 ! print *,"IONS",ions,card(1:6)
4117 else if ((ions).and.(card(1:6).eq.'HETATM')) then
4118 if (firstion.eq.0) then
4122 dc(j,ires)=sccor(j,iii)
4125 call sccenter(ires,iii,sccor)
4128 read (card(12:16),*) atom
4129 print *,"HETATOM", atom(1:2)
4130 read (card(18:20),'(a3)') res
4131 if (atom(3:3).eq.'H') cycle
4132 if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
4133 (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG') &
4134 .or.(atom(1:2).eq.'K ').or.(atom(1:2).eq.'ZN').or.&
4135 (atom(1:2).eq.'O ')) then
4137 print *,"I have water"
4138 if (molecule.ne.5) molecprev=molecule
4140 nres_molec(molecule)=nres_molec(molecule)+1
4141 print *,"HERE",nres_molec(molecule)
4142 if (res.ne.'WAT') res=res(2:3)//' '
4143 itype(ires,molecule)=rescode(ires,res,0,molecule)
4144 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
4148 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
4149 if (ires.eq.0) return
4150 ! Calculate dummy residue coordinates inside the "chain" of a multichain
4153 if (((ires_old.ne.ires).and.(molecule.ne.5)) &
4155 nres_molec(molecule)=nres_molec(molecule)-2
4156 print *,'I have',nres, nres_molec(:)
4158 do k=1,4 ! ions are without dummy
4159 if (nres_molec(k).eq.0) cycle
4161 ! write (iout,*) i,itype(i,1)
4162 ! if (itype(i,1).eq.ntyp1) then
4163 ! write (iout,*) "dummy",i,itype(i,1)
4165 ! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
4166 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2
4170 if (itype(i,k).eq.ntyp1_molec(k)) then
4171 if (itype(i+1,k).eq.ntyp1_molec(k)) then
4172 if (itype(i+2,k).eq.0) then
4173 ! print *,"masz sieczke"
4175 if (itype(i+2,j).ne.0) then
4177 itype(i+1,j)=ntyp1_molec(j)
4178 nres_molec(k)=nres_molec(k)-1
4179 nres_molec(j)=nres_molec(j)+1
4185 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
4186 ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
4187 ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
4188 ! if (unres_pdb) then
4189 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
4190 ! print *,i,'tu dochodze'
4191 ! call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
4199 ! c(j,i)=c(j,i-1)-1.9d0*e2(j)
4203 dcj=(c(j,i-2)-c(j,i-3))/2.0
4204 if (dcj.eq.0) dcj=1.23591524223
4209 else !itype(i+1,1).eq.ntyp1
4210 ! if (unres_pdb) then
4211 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
4212 ! call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
4219 ! c(j,i)=c(j,i+1)-1.9d0*e2(j)
4220 c(j,i)=c(j,i-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
4224 dcj=(c(j,i+3)-c(j,i+2))/2.0
4225 if (dcj.eq.0) dcj=1.23591524223
4230 endif !itype(i+1,1).eq.ntyp1
4231 endif !itype.eq.ntyp1
4235 ! Calculate the CM of the last side chain.
4239 dc(j,ires)=sccor(j,iii)
4242 call sccenter(ires,iii,sccor)
4248 ! print *,"molecule",molecule
4249 if ((itype(nres,1).ne.10)) then
4252 if (molecule.eq.5) molecule=molecprev
4253 itype(nres,molecule)=ntyp1_molec(molecule)
4254 nres_molec(molecule)=nres_molec(molecule)+1
4255 ! if (unres_pdb) then
4256 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
4257 ! call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
4264 ! c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
4268 dcj=(c(j,nres-nres_molec(5)-2)-c(j,nres-nres_molec(5)-3))/2.0
4269 c(j,nres)=c(j,nres-nres_molec(5)-1)+dcj
4270 c(j,2*nres)=c(j,nres)
4274 ! print *,'I have',nres, nres_molec(:)
4276 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
4278 if (nres.ne.nres0) then
4279 write (iout,*) "Error: wrong parameter value: NRES=",nres,&
4281 stop "Error nres value in WHAM input"
4284 !---------------------------------
4285 !el reallocate tables
4288 ! hfrag_alloc(j,i)=hfrag(j,i)
4291 ! bfrag_alloc(j,i)=bfrag(j,i)
4297 ! allocate(hfrag(2,nres/3)) !(2,maxres/3)
4298 !el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
4299 !el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
4300 ! allocate(bfrag(4,nres/3)) !(4,maxres/3)
4304 ! hfrag(j,i)=hfrag_alloc(j,i)
4309 ! bfrag(j,i)=bfrag_alloc(j,i)
4312 !el end reallocate tables
4313 !---------------------------------
4321 c(j,2*nres)=c(j,nres)
4324 if (itype(1,1).eq.ntyp1) then
4328 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
4329 call refsys(2,3,4,e1,e2,e3,fail)
4336 ! c(j,1)=c(j,2)-1.9d0*e2(j)
4337 c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/sqrt(2.0d0)
4341 dcj=(c(j,4)-c(j,3))/2.0
4347 ! First lets assign correct dummy to correct type of chain
4349 if (itype(1,1).eq.ntyp1) then
4350 if (itype(2,1).eq.0) then
4352 if (itype(2,j).ne.0) then
4354 itype(1,j)=ntyp1_molec(j)
4355 nres_molec(1)=nres_molec(1)-1
4356 nres_molec(j)=nres_molec(j)+1
4363 print *,'I have',nres, nres_molec(:)
4365 ! Copy the coordinates to reference coordinates
4371 ! Calculate internal coordinates.
4373 write (iout,'(/a)') &
4374 "Cartesian coordinates of the reference structure"
4375 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
4376 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
4378 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
4379 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
4380 (c(j,ires+nres),j=1,3)
4383 ! znamy już nres wiec mozna alokowac tablice
4384 ! Calculate internal coordinates.
4385 if(me.eq.king.or..not.out1file)then
4386 write (iout,'(a)') &
4387 "Backbone and SC coordinates as read from the PDB"
4389 write (iout,'(i5,i3,2x,a,3f8.3,5x,3f8.3)') &
4390 ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
4391 (c(j,nres+ires),j=1,3)
4394 ! NOW LETS ROCK! SORTING
4395 allocate(c_temporary(3,2*nres))
4396 allocate(itype_temporary(nres,5))
4397 if (.not.allocated(molnum)) allocate(molnum(nres+1))
4398 if (.not.allocated(istype)) write(iout,*) &
4399 "SOMETHING WRONG WITH ISTYTPE"
4400 allocate(istype_temp(nres))
4401 itype_temporary(:,:)=0
4405 if (itype(i,k).ne.0) then
4407 c_temporary(j,seqalingbegin)=c(j,i)
4408 c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
4410 itype_temporary(seqalingbegin,k)=itype(i,k)
4411 print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin
4412 istype_temp(seqalingbegin)=istype(i)
4413 molnum(seqalingbegin)=k
4414 seqalingbegin=seqalingbegin+1
4420 c(j,i)=c_temporary(j,i)
4422 if ((molnum(i-nres).eq.4)) c(j,i)=0.0d0
4428 itype(i,k)=itype_temporary(i,k)
4429 istype(i)=istype_temp(i)
4432 if ((itype(1,1).eq.ntyp1).and.itype(2,5).ne.0) then
4433 ! I have only ions now dummy atoms in the system
4435 itype(1,5)=itype(2,5)
4441 itype(i,5)=itype(i+1,5)
4448 nres_molec(1)=nres_molec(1)-1
4450 ! if (itype(1,1).eq.ntyp1) then
4453 ! if (unres_pdb) then
4454 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
4455 ! call refsys(2,3,4,e1,e2,e3,fail)
4462 ! c(j,1)=c(j,2)-1.9d0*e2(j)
4466 ! dcj=(c(j,4)-c(j,3))/2.0
4468 ! c(j,nres+1)=c(j,1)
4474 write (iout,'(/a)') &
4475 "Cartesian coordinates of the reference structure after sorting"
4476 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
4477 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
4479 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
4480 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
4481 (c(j,ires+nres),j=1,3)
4485 print *,seqalingbegin,nres
4486 if(.not.allocated(vbld)) then
4487 allocate(vbld(2*nres))
4492 if(.not.allocated(vbld_inv)) then
4493 allocate(vbld_inv(2*nres))
4499 if(.not.allocated(theta)) then
4500 allocate(theta(nres+2))
4504 if(.not.allocated(phi)) allocate(phi(nres+2))
4505 if(.not.allocated(alph)) allocate(alph(nres+2))
4506 if(.not.allocated(omeg)) allocate(omeg(nres+2))
4507 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
4508 if(.not.allocated(phiref)) allocate(phiref(nres+2))
4509 if(.not.allocated(costtab)) allocate(costtab(nres))
4510 if(.not.allocated(sinttab)) allocate(sinttab(nres))
4511 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
4512 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
4513 if(.not.allocated(xxref)) allocate(xxref(nres))
4514 if(.not.allocated(yyref)) allocate(yyref(nres))
4515 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
4516 if(.not.allocated(dc_norm)) then
4517 ! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
4518 allocate(dc_norm(3,0:2*nres+2))
4521 write(iout,*) "before int_from_cart"
4522 call int_from_cart(.true.,.false.)
4523 call sc_loc_geom(.false.)
4524 write(iout,*) "after int_from_cart"
4528 thetaref(i)=theta(i)
4531 write(iout,*) "after thetaref"
4539 dc(j,i)=c(j,i+1)-c(j,i)
4540 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
4545 dc(j,i+nres)=c(j,i+nres)-c(j,i)
4546 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
4548 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
4552 ! Copy the coordinates to reference coordinates
4553 ! Splits to single chain if occurs
4557 ! cref(j,i,cou)=c(j,i)
4563 ! chomo(j,i,k)=c(j,i)
4566 ! write(iout,*) "after chomo"
4568 if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
4569 if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
4570 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
4571 !-----------------------------
4575 write (iout,*) "symetr", symetr
4578 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
4580 ! if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
4581 ! chain_length=lll-1
4583 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
4588 cref(j,i,cou)=c(j,i)
4589 cref(j,i+nres,cou)=c(j,i+nres)
4591 chain_rep(j,lll,kkk)=c(j,i)
4592 chain_rep(j,lll+nres,kkk)=c(j,i+nres)
4597 ! write (iout,*) "spraw lancuchy",chain_length,symetr
4599 ! do kkk=1,chain_length
4600 ! write (iout,*) itype(kkk,1),(chain_rep(j,kkk,i), j=1,3)
4604 ! makes copy of chains
4605 write (iout,*) "symetr", symetr
4610 ! if (symetr.gt.1) then
4611 ! call permut(symetr)
4617 ! write(iout,*) (tabperm(i,kkk),kkk=1,4)
4622 ! icha=tabperm(i,kkk)
4623 ! write (iout,*) i,icha
4624 ! do lll=1,chain_length
4626 ! if (cou.le.nres) then
4628 ! kupa=mod(lll,chain_length)
4629 ! iprzes=(kkk-1)*chain_length+lll
4630 ! if (kupa.eq.0) kupa=chain_length
4631 ! write (iout,*) "kupa", kupa
4632 ! cref(j,iprzes,i)=chain_rep(j,kupa,icha)
4633 ! cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
4640 !-koniec robienia kopii
4643 write (iout,*) "nowa struktura", nperm
4645 write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&
4647 cref(3,i,kkk),cref(1,nres+i,kkk),&
4648 cref(2,nres+i,kkk),cref(3,nres+i,kkk)
4650 100 format (//' alpha-carbon coordinates ',&
4651 ' centroid coordinates'/ &
4652 ' ', 6X,'X',11X,'Y',11X,'Z', &
4653 10X,'X',11X,'Y',11X,'Z')
4654 110 format (a,'(',i5,')',6f12.5)
4660 bfrag(i,j)=bfrag(i,j)-ishift
4666 hfrag(i,j)=hfrag(i,j)-ishift
4672 end subroutine readpdb
4673 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
4674 !-----------------------------------------------------------------------------
4676 !-----------------------------------------------------------------------------
4677 subroutine read_control
4691 use random, only: random_init
4692 ! implicit real*8 (a-h,o-z)
4693 ! include 'DIMENSIONS'
4695 use prng, only:prng_restart
4697 logical :: OKRandom!, prng_restart
4700 ! include 'COMMON.IOUNITS'
4701 ! include 'COMMON.TIME1'
4702 ! include 'COMMON.THREAD'
4703 ! include 'COMMON.SBRIDGE'
4704 ! include 'COMMON.CONTROL'
4705 ! include 'COMMON.MCM'
4706 ! include 'COMMON.MAP'
4707 ! include 'COMMON.HEADER'
4708 ! include 'COMMON.CSA'
4709 ! include 'COMMON.CHAIN'
4710 ! include 'COMMON.MUCA'
4711 ! include 'COMMON.MD'
4712 ! include 'COMMON.FFIELD'
4713 ! include 'COMMON.INTERACT'
4714 ! include 'COMMON.SETUP'
4715 !el integer :: KDIAG,ICORFL,IXDR
4716 !el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
4717 character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
4718 'EVVRSP ','Givens ','Jacobi '/),shape(diagmeth))
4719 ! character(len=80) :: ucase
4720 character(len=640) :: controlcard
4722 real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
4724 usampl=.false. ! the default value of usample should be 0
4725 ! write(iout,*) "KURWA2",usampl
4729 read (INP,'(a)') titel
4730 call card_concat(controlcard,.true.)
4731 ! out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
4732 ! print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
4733 call reada(controlcard,'SEED',seed,0.0D0)
4734 call random_init(seed)
4735 ! Set up the time limit (caution! The time must be input in minutes!)
4736 read_cart=index(controlcard,'READ_CART').gt.0
4737 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
4738 call readi(controlcard,'SYM',symetr,1)
4739 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
4740 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
4741 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
4742 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
4743 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
4744 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
4745 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
4746 call reada(controlcard,'DRMS',drms,0.1D0)
4747 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
4748 read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
4749 out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0
4750 out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0
4751 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4752 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
4753 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
4754 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
4755 write (iout,'(a,f10.1)')'DRMS = ',drms
4756 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
4757 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
4759 call readi(controlcard,'NZ_START',nz_start,0)
4760 call readi(controlcard,'NZ_END',nz_end,0)
4761 call readi(controlcard,'IZ_SC',iz_sc,0)
4762 timlim=60.0D0*timlim
4763 safety = 60.0d0*safety
4766 call reada(controlcard,"T_BATH",t_bath,300.0d0)
4767 !C SHIELD keyword sets if the shielding effect of side-chains is used
4768 !C 0 denots no shielding is used all peptide are equally despite the
4769 !C solvent accesible area
4770 !C 1 the newly introduced function
4771 !C 2 reseved for further possible developement
4772 call readi(controlcard,'SHIELD',shield_mode,0)
4773 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4774 write(iout,*) "shield_mode",shield_mode
4775 call readi(controlcard,'TORMODE',tor_mode,0)
4776 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4777 write(iout,*) "torsional and valence angle mode",tor_mode
4779 !C Varibles set size of box
4780 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
4781 protein=index(controlcard,"PROTEIN").gt.0
4782 ions=index(controlcard,"IONS").gt.0
4783 fodson=index(controlcard,"FODSON").gt.0
4784 call readi(controlcard,'OLDION',oldion,1)
4785 nucleic=index(controlcard,"NUCLEIC").gt.0
4786 write (iout,*) "with_theta_constr ",with_theta_constr
4787 AFMlog=(index(controlcard,'AFM'))
4788 selfguide=(index(controlcard,'SELFGUIDE'))
4789 print *,'AFMlog',AFMlog,selfguide,"KUPA"
4790 call readi(controlcard,'GENCONSTR',genconstr,0)
4791 call reada(controlcard,'BOXX',boxxsize,100.0d0)
4792 call reada(controlcard,'BOXY',boxysize,100.0d0)
4793 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
4794 call readi(controlcard,'TUBEMOD',tubemode,0)
4795 print *,"SCELE",scelemode
4796 call readi(controlcard,"SCELEMODE",scelemode,0)
4797 print *,"SCELE",scelemode
4799 ! elemode = 0 is orignal UNRES electrostatics
4800 ! elemode = 1 is "Momo" potentials in progress
4801 ! elemode = 2 is in development EVALD
4804 write (iout,*) TUBEmode,"TUBEMODE"
4805 if (TUBEmode.gt.0) then
4806 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
4807 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
4808 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
4809 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
4810 call reada(controlcard,"VNANO",velnanoconst,0.0d0)
4811 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
4812 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
4813 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
4814 buftubebot=bordtubebot+tubebufthick
4815 buftubetop=bordtubetop-tubebufthick
4818 ! CUTOFFF ON ELECTROSTATICS
4819 call reada(controlcard,"R_CUT_ELE",r_cut_ele,25.0d0)
4820 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
4821 write(iout,*) "R_CUT_ELE=",r_cut_ele
4822 call reada(controlcard,"R_CUT_MART",r_cut_mart,15.0d0)
4823 call reada(controlcard,"LAMBDA_MART",rlamb_mart,0.3d0)
4824 call reada(controlcard,"R_CUT_ANG",r_cut_ang,4.2d0)
4826 ! Lipidec parameters
4827 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
4828 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
4829 if (lipthick.gt.0.0d0) then
4830 bordliptop=(boxzsize+lipthick)/2.0
4831 bordlipbot=bordliptop-lipthick
4832 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
4833 write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
4834 buflipbot=bordlipbot+lipbufthick
4835 bufliptop=bordliptop-lipbufthick
4836 if ((lipbufthick*2.0d0).gt.lipthick) &
4837 write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
4838 endif !lipthick.gt.0
4839 write(iout,*) "bordliptop=",bordliptop
4840 write(iout,*) "bordlipbot=",bordlipbot
4841 write(iout,*) "bufliptop=",bufliptop
4842 write(iout,*) "buflipbot=",buflipbot
4843 write (iout,*) "SHIELD MODE",shield_mode
4845 !C-------------------------
4846 minim=(index(controlcard,'MINIMIZE').gt.0)
4847 dccart=(index(controlcard,'CART').gt.0)
4848 overlapsc=(index(controlcard,'OVERLAP').gt.0)
4849 overlapsc=.not.overlapsc
4850 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
4851 searchsc=.not.searchsc
4852 sideadd=(index(controlcard,'SIDEADD').gt.0)
4853 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
4854 outpdb=(index(controlcard,'PDBOUT').gt.0)
4855 outmol2=(index(controlcard,'MOL2OUT').gt.0)
4856 pdbref=(index(controlcard,'PDBREF').gt.0)
4857 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
4858 indpdb=index(controlcard,'PDBSTART')
4859 extconf=(index(controlcard,'EXTCONF').gt.0)
4860 call readi(controlcard,'IPRINT',iprint,0)
4861 call readi(controlcard,'MAXGEN',maxgen,10000)
4862 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
4863 call readi(controlcard,"KDIAG",kdiag,0)
4864 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
4865 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
4866 write (iout,*) "RESCALE_MODE",rescale_mode
4867 split_ene=index(controlcard,'SPLIT_ENE').gt.0
4868 if (index(controlcard,'REGULAR').gt.0.0D0) then
4869 call reada(controlcard,'WEIDIS',weidis,0.1D0)
4873 if (index(controlcard,'CHECKGRAD').gt.0) then
4875 if (index(controlcard,'CART').gt.0) then
4877 elseif (index(controlcard,'CARINT').gt.0) then
4882 elseif (index(controlcard,'THREAD').gt.0) then
4884 call readi(controlcard,'THREAD',nthread,0)
4885 if (nthread.gt.0) then
4886 call reada(controlcard,'WEIDIS',weidis,0.1D0)
4889 write (iout,'(a)')'A number has to follow the THREAD keyword.'
4890 stop 'Error termination in Read_Control.'
4892 else if (index(controlcard,'MCMA').gt.0) then
4894 else if (index(controlcard,'MCEE').gt.0) then
4896 else if (index(controlcard,'MULTCONF').gt.0) then
4898 else if (index(controlcard,'MAP').gt.0) then
4900 call readi(controlcard,'MAP',nmap,0)
4901 else if (index(controlcard,'CSA').gt.0) then
4903 !rc else if (index(controlcard,'ZSCORE').gt.0) then
4905 !rc ZSCORE is rm from UNRES, modecalc=9 is available
4908 !fcm else if (index(controlcard,'MCMF').gt.0) then
4910 else if (index(controlcard,'SOFTREG').gt.0) then
4912 else if (index(controlcard,'CHECK_BOND').gt.0) then
4914 else if (index(controlcard,'TEST').gt.0) then
4916 else if (index(controlcard,'MD').gt.0) then
4918 else if (index(controlcard,'RE ').gt.0) then
4922 lmuca=index(controlcard,'MUCA').gt.0
4923 call readi(controlcard,'MUCADYN',mucadyn,0)
4924 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
4925 if (lmuca .and. (me.eq.king .or. .not.out1file )) &
4927 write (iout,*) 'MUCADYN=',mucadyn
4928 write (iout,*) 'MUCASMOOTH=',muca_smooth
4931 iscode=index(controlcard,'ONE_LETTER')
4932 indphi=index(controlcard,'PHI')
4933 indback=index(controlcard,'BACK')
4934 iranconf=index(controlcard,'RAND_CONF')
4935 i2ndstr=index(controlcard,'USE_SEC_PRED')
4936 gradout=index(controlcard,'GRADOUT').gt.0
4937 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
4938 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
4939 if (me.eq.king .or. .not.out1file ) &
4940 write (iout,*) "DISTCHAINMAX",distchainmax
4942 if(me.eq.king.or..not.out1file) &
4943 write (iout,'(2a)') diagmeth(kdiag),&
4944 ' routine used to diagonalize matrices.'
4945 if (shield_mode.gt.0) then
4946 pi=4.0D0*datan(1.0D0)
4947 !C VSolvSphere the volume of solving sphere
4949 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
4950 !C there will be no distinction between proline peptide group and normal peptide
4951 !C group in case of shielding parameters
4952 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
4953 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
4954 write (iout,*) VSolvSphere,VSolvSphere_div
4955 !C long axis of side chain
4957 ! long_r_sidechain(i)=vbldsc0(1,i)
4958 ! short_r_sidechain(i)=sigma0(i)
4959 ! write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
4966 end subroutine read_control
4967 !-----------------------------------------------------------------------------
4968 subroutine read_REMDpar
4970 ! Read REMD settings
4977 use control_data, only:out1file
4978 ! implicit real*8 (a-h,o-z)
4979 ! include 'DIMENSIONS'
4980 ! include 'COMMON.IOUNITS'
4981 ! include 'COMMON.TIME1'
4982 ! include 'COMMON.MD'
4985 !el include 'COMMON.LANGEVIN'
4987 !el include 'COMMON.LANGEVIN.lang0'
4989 ! include 'COMMON.INTERACT'
4990 ! include 'COMMON.NAMES'
4991 ! include 'COMMON.GEO'
4992 ! include 'COMMON.REMD'
4993 ! include 'COMMON.CONTROL'
4994 ! include 'COMMON.SETUP'
4995 ! character(len=80) :: ucase
4996 character(len=320) :: controlcard
4997 character(len=3200) :: controlcard1
4998 integer :: iremd_m_total
5001 ! real(kind=8) :: var,ene
5003 if(me.eq.king.or..not.out1file) &
5004 write (iout,*) "REMD setup"
5006 call card_concat(controlcard,.true.)
5007 call readi(controlcard,"NREP",nrep,3)
5008 call readi(controlcard,"NSTEX",nstex,1000)
5009 call reada(controlcard,"RETMIN",retmin,10.0d0)
5010 call reada(controlcard,"RETMAX",retmax,1000.0d0)
5011 mremdsync=(index(controlcard,'SYNC').gt.0)
5012 call readi(controlcard,"NSYN",i_sync_step,100)
5013 restart1file=(index(controlcard,'REST1FILE').gt.0)
5014 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
5015 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
5016 if(max_cache_traj_use.gt.max_cache_traj) &
5017 max_cache_traj_use=max_cache_traj
5018 if(me.eq.king.or..not.out1file) then
5019 !d if (traj1file) then
5020 !rc caching is in testing - NTWX is not ignored
5021 !d write (iout,*) "NTWX value is ignored"
5022 !d write (iout,*) " trajectory is stored to one file by master"
5023 !d write (iout,*) " before exchange at NSTEX intervals"
5025 write (iout,*) "NREP= ",nrep
5026 write (iout,*) "NSTEX= ",nstex
5027 write (iout,*) "SYNC= ",mremdsync
5028 write (iout,*) "NSYN= ",i_sync_step
5029 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
5032 allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
5033 if (index(controlcard,'TLIST').gt.0) then
5035 call card_concat(controlcard1,.true.)
5036 read(controlcard1,*) (remd_t(i),i=1,nrep)
5037 if(me.eq.king.or..not.out1file) &
5038 write (iout,*)'tlist',(remd_t(i),i=1,nrep)
5041 if (index(controlcard,'MLIST').gt.0) then
5043 call card_concat(controlcard1,.true.)
5044 read(controlcard1,*) (remd_m(i),i=1,nrep)
5045 if(me.eq.king.or..not.out1file) then
5046 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
5049 iremd_m_total=iremd_m_total+remd_m(i)
5051 write (iout,*) 'Total number of replicas ',iremd_m_total
5054 if(me.eq.king.or..not.out1file) &
5055 write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
5057 end subroutine read_REMDpar
5058 !-----------------------------------------------------------------------------
5059 subroutine read_MDpar
5063 use control_data, only: r_cut,rlamb,out1file,r_cut_mart,rlamb_mart
5065 use geometry_data, only: pi
5067 ! implicit real*8 (a-h,o-z)
5068 ! include 'DIMENSIONS'
5069 ! include 'COMMON.IOUNITS'
5070 ! include 'COMMON.TIME1'
5071 ! include 'COMMON.MD'
5074 !el include 'COMMON.LANGEVIN'
5076 !el include 'COMMON.LANGEVIN.lang0'
5078 ! include 'COMMON.INTERACT'
5079 ! include 'COMMON.NAMES'
5080 ! include 'COMMON.GEO'
5081 ! include 'COMMON.SETUP'
5082 ! include 'COMMON.CONTROL'
5083 ! include 'COMMON.SPLITELE'
5084 ! character(len=80) :: ucase
5085 character(len=320) :: controlcard
5090 call card_concat(controlcard,.true.)
5091 call readi(controlcard,"NSTEP",n_timestep,1000000)
5092 call readi(controlcard,"NTWE",ntwe,100)
5093 call readi(controlcard,"NTWX",ntwx,1000)
5094 call readi(controlcard,"NFOD",nfodstep,100)
5095 call reada(controlcard,"DT",d_time,1.0d-1)
5096 call reada(controlcard,"DVMAX",dvmax,2.0d1)
5097 call reada(controlcard,"DAMAX",damax,1.0d1)
5098 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
5099 call readi(controlcard,"LANG",lang,0)
5100 RESPA = index(controlcard,"RESPA") .gt. 0
5101 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
5102 ntime_split0=ntime_split
5103 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
5104 ntime_split0=ntime_split
5105 call reada(controlcard,"R_CUT",r_cut,2.0d0)
5106 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
5107 rest = index(controlcard,"REST").gt.0
5108 tbf = index(controlcard,"TBF").gt.0
5109 usampl = index(controlcard,"USAMPL").gt.0
5110 ! write(iout,*) "KURWA",usampl
5111 mdpdb = index(controlcard,"MDPDB").gt.0
5112 call reada(controlcard,"T_BATH",t_bath,300.0d0)
5113 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
5114 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
5115 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
5116 if (count_reset_moment.eq.0) count_reset_moment=1000000000
5117 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
5118 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
5119 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
5120 if (count_reset_vel.eq.0) count_reset_vel=1000000000
5121 large = index(controlcard,"LARGE").gt.0
5122 print_compon = index(controlcard,"PRINT_COMPON").gt.0
5123 rattle = index(controlcard,"RATTLE").gt.0
5124 preminim=(index(controlcard,'PREMINIM').gt.0)
5125 forceminim=(index(controlcard,'FORCEMINIM').gt.0)
5126 write (iout,*) "PREMINIM ",preminim
5127 dccart=(index(controlcard,'CART').gt.0)
5128 if (preminim) call read_minim
5129 ! if performing umbrella sampling, fragments constrained are read from the fragment file
5135 if(me.eq.king.or..not.out1file) then
5137 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
5139 write (iout,'(a)') "The units are:"
5140 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
5141 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
5142 " acceleration: angstrom/(48.9 fs)**2"
5143 write (iout,'(a)') "energy: kcal/mol, temperature: K"
5145 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
5146 write (iout,'(a60,f10.5,a)') &
5147 "Initial time step of numerical integration:",d_time,&
5149 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
5151 write (iout,'(2a,i4,a)') &
5152 "A-MTS algorithm used; initial time step for fast-varying",&
5153 " short-range forces split into",ntime_split," steps."
5154 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
5155 r_cut," lambda",rlamb
5157 write (iout,'(2a,f10.5)') &
5158 "Maximum acceleration threshold to reduce the time step",&
5159 "/increase split number:",damax
5160 write (iout,'(2a,f10.5)') &
5161 "Maximum predicted energy drift to reduce the timestep",&
5162 "/increase split number:",edriftmax
5163 write (iout,'(a60,f10.5)') &
5164 "Maximum velocity threshold to reduce velocities:",dvmax
5165 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
5166 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
5167 if (rattle) write (iout,'(a60)') &
5168 "Rattle algorithm used to constrain the virtual bonds"
5172 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
5173 call reada(controlcard,"RWAT",rwat,1.4d0)
5174 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
5175 surfarea=index(controlcard,"SURFAREA").gt.0
5176 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
5177 if(me.eq.king.or..not.out1file)then
5178 write (iout,'(/a,$)') "Langevin dynamics calculation"
5180 write (iout,'(a/)') &
5181 " with direct integration of Langevin equations"
5182 else if (lang.eq.2) then
5183 write (iout,'(a/)') " with TINKER stochasic MD integrator"
5184 else if (lang.eq.3) then
5185 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
5186 else if (lang.eq.4) then
5187 write (iout,'(a/)') " in overdamped mode"
5189 write (iout,'(//a,i5)') &
5190 "=========== ERROR: Unknown Langevin dynamics mode:",lang
5193 write (iout,'(a60,f10.5)') "Temperature:",t_bath
5194 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
5195 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
5196 write (iout,'(a60,f10.5)') &
5197 "Scaling factor of the friction forces:",scal_fric
5198 if (surfarea) write (iout,'(2a,i10,a)') &
5199 "Friction coefficients will be scaled by solvent-accessible",&
5200 " surface area every",reset_fricmat," steps."
5202 ! Calculate friction coefficients and bounds of stochastic forces
5203 eta=6*pi*cPoise*etawat
5204 if(me.eq.king.or..not.out1file) &
5205 write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
5208 do j=1,5 !types of molecules
5209 gamp(j)=scal_fric*(pstok(j)+rwat)*eta
5210 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
5212 allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1)
5213 do j=1,5 !types of molecules
5215 gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta
5216 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
5220 if(me.eq.king.or..not.out1file)then
5222 write (iout,'(/2a/)') &
5223 "Radii of site types and friction coefficients and std's of",&
5224 " stochastic forces of fully exposed sites"
5225 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(j),stdfp*dsqrt(gamp(j))
5228 write (iout,'(a5,f5.2,2f10.5)') restyp(i,j),restok(i,j),&
5229 gamsc(i,j),stdfsc(i,j)*dsqrt(gamsc(i,j))
5234 if(me.eq.king.or..not.out1file)then
5235 write (iout,'(a)') "Berendsen bath calculation"
5236 write (iout,'(a60,f10.5)') "Temperature:",t_bath
5237 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
5239 write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
5240 count_reset_moment," steps"
5242 write (iout,'(a,i10,a)') &
5243 "Velocities will be reset at random every",count_reset_vel,&
5247 if(me.eq.king.or..not.out1file) &
5248 write (iout,'(a31)') "Microcanonical mode calculation"
5250 if(me.eq.king.or..not.out1file)then
5251 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
5253 write(iout,*) "MD running with constraints."
5254 write(iout,*) "Equilibration time ", eq_time, " mtus."
5255 write(iout,*) "Constraining ", nfrag," fragments."
5256 write(iout,*) "Length of each fragment, weight and q0:"
5258 write (iout,*) "Set of restraints #",iset
5260 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
5261 ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
5263 write(iout,*) "constraints between ", npair, "fragments."
5264 write(iout,*) "constraint pairs, weights and q0:"
5266 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
5267 ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
5269 write(iout,*) "angle constraints within ", nfrag_back,&
5270 "backbone fragments."
5271 write(iout,*) "fragment, weights:"
5273 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
5274 ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
5275 wfrag_back(2,i,iset),wfrag_back(3,i,iset)
5278 iset=mod(kolor,nset)+1
5281 if(me.eq.king.or..not.out1file) &
5282 write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
5284 end subroutine read_MDpar
5285 !-----------------------------------------------------------------------------
5289 ! implicit real*8 (a-h,o-z)
5290 ! include 'DIMENSIONS'
5291 ! include 'COMMON.MAP'
5292 ! include 'COMMON.IOUNITS'
5293 character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
5294 character(len=80) :: mapcard !,ucase
5297 ! real(kind=8) :: var,ene
5300 read (inp,'(a)') mapcard
5301 mapcard=ucase(mapcard)
5302 if (index(mapcard,'PHI').gt.0) then
5304 else if (index(mapcard,'THE').gt.0) then
5306 else if (index(mapcard,'ALP').gt.0) then
5308 else if (index(mapcard,'OME').gt.0) then
5311 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
5312 stop 'Error - illegal variable spec in MAP card.'
5314 call readi (mapcard,'RES1',res1(imap),0)
5315 call readi (mapcard,'RES2',res2(imap),0)
5316 if (res1(imap).eq.0) then
5317 res1(imap)=res2(imap)
5318 else if (res2(imap).eq.0) then
5319 res2(imap)=res1(imap)
5321 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
5322 write (iout,'(a)') &
5323 'Error - illegal definition of variable group in MAP.'
5324 stop 'Error - illegal definition of variable group in MAP.'
5326 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
5327 call reada(mapcard,'TO',ang_to(imap),0.0D0)
5328 call readi(mapcard,'NSTEP',nstep(imap),0)
5329 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
5330 write (iout,'(a)') &
5331 'Illegal boundary and/or step size specification in MAP.'
5332 stop 'Illegal boundary and/or step size specification in MAP.'
5336 end subroutine map_read
5337 !-----------------------------------------------------------------------------
5340 use control_data, only: vdisulf
5342 ! implicit real*8 (a-h,o-z)
5343 ! include 'DIMENSIONS'
5344 ! include 'COMMON.IOUNITS'
5345 ! include 'COMMON.GEO'
5346 ! include 'COMMON.CSA'
5347 ! include 'COMMON.BANK'
5348 ! include 'COMMON.CONTROL'
5349 ! character(len=80) :: ucase
5350 character(len=620) :: mcmcard
5352 ! integer :: ntf,ik,iw_pdb
5353 ! real(kind=8) :: var,ene
5355 call card_concat(mcmcard,.true.)
5357 call readi(mcmcard,'NCONF',nconf,50)
5358 call readi(mcmcard,'NADD',nadd,0)
5359 call readi(mcmcard,'JSTART',jstart,1)
5360 call readi(mcmcard,'JEND',jend,1)
5361 call readi(mcmcard,'NSTMAX',nstmax,500000)
5362 call readi(mcmcard,'N0',n0,1)
5363 call readi(mcmcard,'N1',n1,6)
5364 call readi(mcmcard,'N2',n2,4)
5365 call readi(mcmcard,'N3',n3,0)
5366 call readi(mcmcard,'N4',n4,0)
5367 call readi(mcmcard,'N5',n5,0)
5368 call readi(mcmcard,'N6',n6,10)
5369 call readi(mcmcard,'N7',n7,0)
5370 call readi(mcmcard,'N8',n8,0)
5371 call readi(mcmcard,'N9',n9,0)
5372 call readi(mcmcard,'N14',n14,0)
5373 call readi(mcmcard,'N15',n15,0)
5374 call readi(mcmcard,'N16',n16,0)
5375 call readi(mcmcard,'N17',n17,0)
5376 call readi(mcmcard,'N18',n18,0)
5378 vdisulf=(index(mcmcard,'DYNSS').gt.0)
5380 call readi(mcmcard,'NDIFF',ndiff,2)
5381 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
5382 call readi(mcmcard,'IS1',is1,1)
5383 call readi(mcmcard,'IS2',is2,8)
5384 call readi(mcmcard,'NRAN0',nran0,4)
5385 call readi(mcmcard,'NRAN1',nran1,2)
5386 call readi(mcmcard,'IRR',irr,1)
5387 call readi(mcmcard,'NSEED',nseed,20)
5388 call readi(mcmcard,'NTOTAL',ntotal,10000)
5389 call reada(mcmcard,'CUT1',cut1,2.0d0)
5390 call reada(mcmcard,'CUT2',cut2,5.0d0)
5391 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
5392 call readi(mcmcard,'ICMAX',icmax,3)
5393 call readi(mcmcard,'IRESTART',irestart,0)
5394 !!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
5397 call reada(mcmcard,'DELE',dele,20.0d0)
5398 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
5399 call readi(mcmcard,'IREF',iref,0)
5400 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
5401 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
5402 call readi(mcmcard,'NCONF_IN',nconf_in,0)
5403 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
5404 write (iout,*) "NCONF_IN",nconf_in
5406 end subroutine csaread
5407 !-----------------------------------------------------------------------------
5411 use control_data, only: MaxMoveType
5414 ! implicit real*8 (a-h,o-z)
5415 ! include 'DIMENSIONS'
5416 ! include 'COMMON.MCM'
5417 ! include 'COMMON.MCE'
5418 ! include 'COMMON.IOUNITS'
5419 ! character(len=80) :: ucase
5420 character(len=320) :: mcmcard
5423 ! real(kind=8) :: var,ene
5425 call card_concat(mcmcard,.true.)
5426 call readi(mcmcard,'MAXACC',maxacc,100)
5427 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
5428 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
5429 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
5430 call readi(mcmcard,'MAXREPM',maxrepm,200)
5431 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
5432 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
5433 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
5434 call reada(mcmcard,'E_UP',e_up,5.0D0)
5435 call reada(mcmcard,'DELTE',delte,0.1D0)
5436 call readi(mcmcard,'NSWEEP',nsweep,5)
5437 call readi(mcmcard,'NSTEPH',nsteph,0)
5438 call readi(mcmcard,'NSTEPC',nstepc,0)
5439 call reada(mcmcard,'TMIN',tmin,298.0D0)
5440 call reada(mcmcard,'TMAX',tmax,298.0D0)
5441 call readi(mcmcard,'NWINDOW',nwindow,0)
5442 call readi(mcmcard,'PRINT_MC',print_mc,0)
5443 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
5444 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
5445 ent_read=(index(mcmcard,'ENT_READ').gt.0)
5446 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
5447 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
5448 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
5449 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
5450 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
5451 if (nwindow.gt.0) then
5452 allocate(winstart(nwindow)) !!el (maxres)
5453 allocate(winend(nwindow)) !!el
5454 allocate(winlen(nwindow)) !!el
5455 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
5457 winlen(i)=winend(i)-winstart(i)+1
5460 if (tmax.lt.tmin) tmax=tmin
5461 if (tmax.eq.tmin) then
5465 if (nstepc.gt.0 .and. nsteph.gt.0) then
5466 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
5467 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
5469 allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
5470 ! Probabilities of different move types
5471 sumpro_type(0)=0.0D0
5472 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
5473 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
5474 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
5475 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
5476 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
5477 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
5478 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
5480 print *,'i',i,' sumprotype',sumpro_type(i)
5481 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
5482 print *,'i',i,' sumprotype',sumpro_type(i)
5485 end subroutine mcmread
5486 !-----------------------------------------------------------------------------
5487 subroutine read_minim
5490 ! implicit real*8 (a-h,o-z)
5491 ! include 'DIMENSIONS'
5492 ! include 'COMMON.MINIM'
5493 ! include 'COMMON.IOUNITS'
5494 ! character(len=80) :: ucase
5495 character(len=320) :: minimcard
5497 ! integer :: ntf,ik,iw_pdb
5498 ! real(kind=8) :: var,ene
5500 call card_concat(minimcard,.true.)
5501 call readi(minimcard,'MAXMIN',maxmin,2000)
5502 call readi(minimcard,'MAXFUN',maxfun,5000)
5503 call readi(minimcard,'MINMIN',minmin,maxmin)
5504 call readi(minimcard,'MINFUN',minfun,maxmin)
5505 call reada(minimcard,'TOLF',tolf,1.0D-2)
5506 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
5507 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
5508 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
5509 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
5510 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
5511 'Options in energy minimization:'
5512 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
5513 'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
5514 'MinMin:',MinMin,' MinFun:',MinFun,&
5515 ' TolF:',TolF,' RTolF:',RTolF
5517 end subroutine read_minim
5518 !-----------------------------------------------------------------------------
5519 subroutine openunits
5521 use MD_data, only: usampl
5524 use control_data, only:out1file
5525 use control, only: getenv_loc
5526 ! implicit real*8 (a-h,o-z)
5527 ! include 'DIMENSIONS'
5530 character(len=16) :: form,nodename
5531 integer :: nodelen,ierror,npos
5533 ! include 'COMMON.SETUP'
5534 ! include 'COMMON.IOUNITS'
5535 ! include 'COMMON.MD'
5536 ! include 'COMMON.CONTROL'
5537 integer :: lenpre,lenpot,lentmp !,ilen
5539 character(len=3) :: out1file_text !,ucase
5540 character(len=3) :: ll
5543 ! integer :: ntf,ik,iw_pdb
5544 ! real(kind=8) :: var,ene
5546 ! print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
5547 call getenv_loc("PREFIX",prefix)
5549 call getenv_loc("POT",pot)
5550 call getenv_loc("DIRTMP",tmpdir)
5551 call getenv_loc("CURDIR",curdir)
5552 call getenv_loc("OUT1FILE",out1file_text)
5553 ! print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
5554 out1file_text=ucase(out1file_text)
5555 if (out1file_text(1:1).eq."Y") then
5558 out1file=fg_rank.gt.0
5563 if (lentmp.gt.0) then
5564 write (*,'(80(1h!))')
5565 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
5566 write (*,'(80(1h!))')
5567 write (*,*)"All output files will be on node /tmp directory."
5569 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
5570 if (me.eq.king) then
5571 write (*,*) "The master node is ",nodename
5572 else if (fg_rank.eq.0) then
5573 write (*,*) "I am the CG slave node ",nodename
5575 write (*,*) "I am the FG slave node ",nodename
5578 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
5579 lenpre = lentmp+lenpre+1
5581 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
5582 ! Get the names and open the input files
5583 #if defined(WINIFL) || defined(WINPGI)
5584 open(1,file=pref_orig(:ilen(pref_orig))// &
5585 '.inp',status='old',readonly,shared)
5586 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5587 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5588 ! Get parameter filenames and open the parameter files.
5589 call getenv_loc('BONDPAR',bondname)
5590 open (ibond,file=bondname,status='old',readonly,shared)
5591 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5592 open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared)
5593 call getenv_loc('THETPAR',thetname)
5594 open (ithep,file=thetname,status='old',readonly,shared)
5595 call getenv_loc('ROTPAR',rotname)
5596 open (irotam,file=rotname,status='old',readonly,shared)
5597 call getenv_loc('TORPAR',torname)
5598 open (itorp,file=torname,status='old',readonly,shared)
5599 call getenv_loc('TORDPAR',tordname)
5600 open (itordp,file=tordname,status='old',readonly,shared)
5601 call getenv_loc('FOURIER',fouriername)
5602 open (ifourier,file=fouriername,status='old',readonly,shared)
5603 call getenv_loc('ELEPAR',elename)
5604 open (ielep,file=elename,status='old',readonly,shared)
5605 call getenv_loc('SIDEPAR',sidename)
5606 open (isidep,file=sidename,status='old',readonly,shared)
5608 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5609 open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared)
5610 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5611 open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared)
5612 call getenv_loc('TORPAR_NUCL',torname_nucl)
5613 open (itorp_nucl,file=torname_nucl,status='old',readonly,shared)
5614 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5615 open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared)
5616 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5617 open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared)
5620 #elif (defined CRAY) || (defined AIX)
5621 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
5623 ! print *,"Processor",myrank," opened file 1"
5624 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5625 ! print *,"Processor",myrank," opened file 9"
5626 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5627 ! Get parameter filenames and open the parameter files.
5628 call getenv_loc('BONDPAR',bondname)
5629 open (ibond,file=bondname,status='old',action='read')
5630 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5631 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
5633 ! print *,"Processor",myrank," opened file IBOND"
5634 call getenv_loc('THETPAR',thetname)
5635 open (ithep,file=thetname,status='old',action='read')
5636 ! print *,"Processor",myrank," opened file ITHEP"
5637 call getenv_loc('ROTPAR',rotname)
5638 open (irotam,file=rotname,status='old',action='read')
5640 call getenv_loc('ROTPAR_END',rotname_end)
5641 open (irotam_end,file=rotname_end,status='old',action='read')
5643 ! print *,"Processor",myrank," opened file IROTAM"
5644 call getenv_loc('TORPAR',torname)
5645 open (itorp,file=torname,status='old',action='read')
5646 ! print *,"Processor",myrank," opened file ITORP"
5647 call getenv_loc('TORDPAR',tordname)
5648 open (itordp,file=tordname,status='old',action='read')
5649 ! print *,"Processor",myrank," opened file ITORDP"
5650 call getenv_loc('SCCORPAR',sccorname)
5651 open (isccor,file=sccorname,status='old',action='read')
5652 ! print *,"Processor",myrank," opened file ISCCOR"
5653 call getenv_loc('FOURIER',fouriername)
5654 open (ifourier,file=fouriername,status='old',action='read')
5655 ! print *,"Processor",myrank," opened file IFOURIER"
5656 call getenv_loc('ELEPAR',elename)
5657 open (ielep,file=elename,status='old',action='read')
5658 ! print *,"Processor",myrank," opened file IELEP"
5659 call getenv_loc('SIDEPAR',sidename)
5660 open (isidep,file=sidename,status='old',action='read')
5662 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5663 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
5664 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5665 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
5666 call getenv_loc('TORPAR_NUCL',torname_nucl)
5667 open (itorp_nucl,file=torname_nucl,status='old',action='read')
5668 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5669 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
5670 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5671 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
5673 call getenv_loc('LIPTRANPAR',liptranname)
5674 open (iliptranpar,file=liptranname,status='old',action='read')
5675 call getenv_loc('TUBEPAR',tubename)
5676 open (itube,file=tubename,status='old',action='read')
5677 call getenv_loc('IONPAR',ionname)
5678 open (iion,file=ionname,status='old',action='read')
5679 call getenv_loc('IONPAR_TRAN',iontranname)
5680 open (iiontran,file=iontranname,status='old',action='read')
5681 ! print *,"Processor",myrank," opened file ISIDEP"
5682 ! print *,"Processor",myrank," opened parameter files"
5684 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
5685 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5686 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5687 ! Get parameter filenames and open the parameter files.
5688 call getenv_loc('BONDPAR',bondname)
5689 open (ibond,file=bondname,status='old')
5690 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5691 open (ibond_nucl,file=bondname_nucl,status='old')
5693 call getenv_loc('THETPAR',thetname)
5694 open (ithep,file=thetname,status='old')
5695 call getenv_loc('ROTPAR',rotname)
5696 open (irotam,file=rotname,status='old')
5698 call getenv_loc('ROTPAR_END',rotname_end)
5699 open (irotam_end,file=rotname_end,status='old')
5701 call getenv_loc('TORPAR',torname)
5702 open (itorp,file=torname,status='old')
5703 call getenv_loc('TORDPAR',tordname)
5704 open (itordp,file=tordname,status='old')
5705 call getenv_loc('SCCORPAR',sccorname)
5706 open (isccor,file=sccorname,status='old')
5707 call getenv_loc('FOURIER',fouriername)
5708 open (ifourier,file=fouriername,status='old')
5709 call getenv_loc('ELEPAR',elename)
5710 open (ielep,file=elename,status='old')
5711 call getenv_loc('SIDEPAR',sidename)
5712 open (isidep,file=sidename,status='old')
5714 open (ithep_nucl,file=thetname_nucl,status='old')
5715 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5716 open (irotam_nucl,file=rotname_nucl,status='old')
5717 call getenv_loc('TORPAR_NUCL',torname_nucl)
5718 open (itorp_nucl,file=torname_nucl,status='old')
5719 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5720 open (itordp_nucl,file=tordname_nucl,status='old')
5721 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5722 open (isidep_nucl,file=sidename_nucl,status='old')
5724 call getenv_loc('LIPTRANPAR',liptranname)
5725 open (iliptranpar,file=liptranname,status='old')
5726 call getenv_loc('TUBEPAR',tubename)
5727 open (itube,file=tubename,status='old')
5728 call getenv_loc('IONPAR',ionname)
5729 open (iion,file=ionname,status='old')
5730 call getenv_loc('IONPAR_NUCL',ionnuclname)
5731 open (iionnucl,file=ionnuclname,status='old')
5732 call getenv_loc('IONPAR_TRAN',iontranname)
5733 open (iiontran,file=iontranname,status='old')
5734 call getenv_loc('WATWAT',iwaterwatername)
5735 open (iwaterwater,file=iwaterwatername,status='old')
5736 call getenv_loc('WATPROT',iwaterscname)
5737 open (iwatersc,file=iwaterscname,status='old')
5740 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
5742 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5743 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5744 ! Get parameter filenames and open the parameter files.
5745 call getenv_loc('BONDPAR',bondname)
5746 open (ibond,file=bondname,status='old',action='read')
5747 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5748 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
5749 call getenv_loc('THETPAR',thetname)
5750 open (ithep,file=thetname,status='old',action='read')
5751 call getenv_loc('ROTPAR',rotname)
5752 open (irotam,file=rotname,status='old',action='read')
5754 call getenv_loc('ROTPAR_END',rotname_end)
5755 open (irotam_end,file=rotname_end,status='old',action='read')
5757 call getenv_loc('TORPAR',torname)
5758 open (itorp,file=torname,status='old',action='read')
5759 call getenv_loc('TORDPAR',tordname)
5760 open (itordp,file=tordname,status='old',action='read')
5761 call getenv_loc('SCCORPAR',sccorname)
5762 open (isccor,file=sccorname,status='old',action='read')
5764 call getenv_loc('THETPARPDB',thetname_pdb)
5765 print *,"thetname_pdb ",thetname_pdb
5766 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
5767 print *,ithep_pdb," opened"
5769 call getenv_loc('FOURIER',fouriername)
5770 open (ifourier,file=fouriername,status='old',readonly)
5771 call getenv_loc('ELEPAR',elename)
5772 open (ielep,file=elename,status='old',readonly)
5773 call getenv_loc('SIDEPAR',sidename)
5774 open (isidep,file=sidename,status='old',readonly)
5776 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5777 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
5778 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5779 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
5780 call getenv_loc('TORPAR_NUCL',torname_nucl)
5781 open (itorp_nucl,file=torname_nucl,status='old',action='read')
5782 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5783 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
5784 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5785 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
5786 call getenv_loc('SIDEPAR_SCBASE',sidename_scbase)
5787 open (isidep_scbase,file=sidename_scbase,status='old',action='read')
5788 call getenv_loc('PEPPAR_PEPBASE',pepname_pepbase)
5789 open (isidep_pepbase,file=pepname_pepbase,status='old',action='read')
5790 call getenv_loc('SCPAR_PHOSPH',pepname_scpho)
5791 open (isidep_scpho,file=pepname_scpho,status='old',action='read')
5792 call getenv_loc('PEPPAR_PHOSPH',pepname_peppho)
5793 open (isidep_peppho,file=pepname_peppho,status='old',action='read')
5796 call getenv_loc('LIPTRANPAR',liptranname)
5797 open (iliptranpar,file=liptranname,status='old',action='read')
5798 call getenv_loc('LIPBOND',lipbondname)
5799 open (ilipbond,file=lipbondname,status='old',action='read')
5800 call getenv_loc('LIPNONBOND',lipnonbondname)
5801 open (ilipnonbond,file=lipnonbondname,status='old',action='read')
5802 call getenv_loc('TUBEPAR',tubename)
5803 open (itube,file=tubename,status='old',action='read')
5804 call getenv_loc('IONPAR',ionname)
5805 open (iion,file=ionname,status='old',action='read')
5806 call getenv_loc('IONPAR_NUCL',ionnuclname)
5807 open (iionnucl,file=ionnuclname,status='old',action='read')
5808 call getenv_loc('IONPAR_TRAN',iontranname)
5809 open (iiontran,file=iontranname,status='old',action='read')
5810 call getenv_loc('WATWAT',iwaterwatername)
5811 open (iwaterwater,file=iwaterwatername,status='old',action='read')
5812 call getenv_loc('WATPROT',iwaterscname)
5813 open (iwatersc,file=iwaterscname,status='old',action='read')
5816 call getenv_loc('ROTPARPDB',rotname_pdb)
5817 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
5820 call getenv_loc('SCPPAR_NUCL',scpname_nucl)
5821 #if defined(WINIFL) || defined(WINPGI)
5822 open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared)
5823 #elif (defined CRAY) || (defined AIX)
5824 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
5826 open (iscpp_nucl,file=scpname_nucl,status='old')
5828 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
5833 ! 8/9/01 In the newest version SCp interaction constants are read from a file
5834 ! Use -DOLDSCP to use hard-coded constants instead.
5836 call getenv_loc('SCPPAR',scpname)
5837 #if defined(WINIFL) || defined(WINPGI)
5838 open (iscpp,file=scpname,status='old',readonly,shared)
5839 #elif (defined CRAY) || (defined AIX)
5840 open (iscpp,file=scpname,status='old',action='read')
5842 open (iscpp,file=scpname,status='old')
5844 open (iscpp,file=scpname,status='old',action='read')
5847 call getenv_loc('PATTERN',patname)
5848 #if defined(WINIFL) || defined(WINPGI)
5849 open (icbase,file=patname,status='old',readonly,shared)
5850 #elif (defined CRAY) || (defined AIX)
5851 open (icbase,file=patname,status='old',action='read')
5853 open (icbase,file=patname,status='old')
5855 open (icbase,file=patname,status='old',action='read')
5858 ! Open output file only for CG processes
5859 ! print *,"Processor",myrank," fg_rank",fg_rank
5860 if (fg_rank.eq.0) then
5862 if (nodes.eq.1) then
5865 npos = dlog10(dfloat(nodes-1))+1
5867 if (npos.lt.3) npos=3
5868 write (liczba,'(i1)') npos
5869 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
5871 write (liczba,form) me
5872 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
5873 liczba(:ilen(liczba))
5874 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
5876 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
5878 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
5879 liczba(:ilen(liczba))//'.mol2'
5880 statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
5881 liczba(:ilen(liczba))//'.stat'
5883 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
5884 //liczba(:ilen(liczba))//'.stat')
5885 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
5888 qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
5889 liczba(:ilen(liczba))//'.const'
5894 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
5895 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
5896 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
5897 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
5898 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
5900 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
5902 rest2name=prefix(:ilen(prefix))//'.rst'
5904 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
5907 #if defined(AIX) || defined(PGI)
5908 if (me.eq.king .or. .not. out1file) &
5909 open(iout,file=outname,status='unknown')
5911 if (fg_rank.gt.0) then
5912 write (liczba,'(i3.3)') myrank/nfgtasks
5913 write (ll,'(bz,i3.3)') fg_rank
5914 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
5919 open(igeom,file=intname,status='unknown',position='append')
5920 open(ipdb,file=pdbname,status='unknown')
5921 open(imol2,file=mol2name,status='unknown')
5922 open(istat,file=statname,status='unknown',position='append')
5924 !1out open(iout,file=outname,status='unknown')
5927 if (me.eq.king .or. .not.out1file) &
5928 open(iout,file=outname,status='unknown')
5930 if (fg_rank.gt.0) then
5931 write (liczba,'(i3.3)') myrank/nfgtasks
5932 write (ll,'(bz,i3.3)') fg_rank
5933 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
5938 open(igeom,file=intname,status='unknown',access='append')
5939 open(ipdb,file=pdbname,status='unknown')
5940 open(imol2,file=mol2name,status='unknown')
5941 open(istat,file=statname,status='unknown',access='append')
5943 !1out open(iout,file=outname,status='unknown')
5946 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
5947 csa_seed=prefix(:lenpre)//'.CSA.seed'
5948 csa_history=prefix(:lenpre)//'.CSA.history'
5949 csa_bank=prefix(:lenpre)//'.CSA.bank'
5950 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
5951 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
5952 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
5953 !!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
5954 csa_int=prefix(:lenpre)//'.int'
5955 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
5956 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
5957 csa_in=prefix(:lenpre)//'.CSA.in'
5958 ! print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
5961 write (iout,'(80(1h-))')
5962 write (iout,'(30x,a)') "FILE ASSIGNMENT"
5963 write (iout,'(80(1h-))')
5964 write (iout,*) "Input file : ",&
5965 pref_orig(:ilen(pref_orig))//'.inp'
5966 write (iout,*) "Output file : ",&
5967 outname(:ilen(outname))
5969 write (iout,*) "Sidechain potential file : ",&
5970 sidename(:ilen(sidename))
5972 write (iout,*) "SCp potential file : ",&
5973 scpname(:ilen(scpname))
5975 write (iout,*) "Electrostatic potential file : ",&
5976 elename(:ilen(elename))
5977 write (iout,*) "Cumulant coefficient file : ",&
5978 fouriername(:ilen(fouriername))
5979 write (iout,*) "Torsional parameter file : ",&
5980 torname(:ilen(torname))
5981 write (iout,*) "Double torsional parameter file : ",&
5982 tordname(:ilen(tordname))
5983 write (iout,*) "SCCOR parameter file : ",&
5984 sccorname(:ilen(sccorname))
5985 write (iout,*) "Bond & inertia constant file : ",&
5986 bondname(:ilen(bondname))
5987 write (iout,*) "Bending parameter file : ",&
5988 thetname(:ilen(thetname))
5989 write (iout,*) "Rotamer parameter file : ",&
5990 rotname(:ilen(rotname))
5993 write (iout,*) "Thetpdb parameter file : ",&
5994 thetname_pdb(:ilen(thetname_pdb))
5997 write (iout,*) "Threading database : ",&
5998 patname(:ilen(patname))
6000 write (iout,*)" DIRTMP : ",&
6002 write (iout,'(80(1h-))')
6005 end subroutine openunits
6006 !-----------------------------------------------------------------------------
6009 use geometry_data, only: nres,dc
6011 ! implicit real*8 (a-h,o-z)
6012 ! include 'DIMENSIONS'
6013 ! include 'COMMON.CHAIN'
6014 ! include 'COMMON.IOUNITS'
6015 ! include 'COMMON.MD'
6018 ! real(kind=8) :: var,ene
6020 open(irest2,file=rest2name,status='unknown')
6021 read(irest2,*) totT,EK,potE,totE,t_bath
6024 ! AL 4/17/17: Now reading d_t(0,:) too
6026 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
6029 ! AL 4/17/17: Now reading d_c(0,:) too
6031 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
6034 read (irest2,*) iset
6038 end subroutine readrst
6039 !-----------------------------------------------------------------------------
6040 subroutine read_fragments
6044 use control_data, only:out1file
6047 ! implicit real*8 (a-h,o-z)
6048 ! include 'DIMENSIONS'
6052 ! include 'COMMON.SETUP'
6053 ! include 'COMMON.CHAIN'
6054 ! include 'COMMON.IOUNITS'
6055 ! include 'COMMON.MD'
6056 ! include 'COMMON.CONTROL'
6059 ! real(kind=8) :: var,ene
6061 read(inp,*) nset,nfrag,npair,nfrag_back
6063 !el from module energy
6064 ! if(.not.allocated(mset)) allocate(mset(nset)) !(maxprocs/20)
6065 if(.not.allocated(wfrag_back)) then
6066 allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
6067 allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
6069 allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
6070 allocate(ifrag(2,nfrag,nset)) !(2,50,maxprocs/20)
6072 allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
6073 allocate(ipair(2,npair,nset)) !(2,100,maxprocs/20)
6076 if(me.eq.king.or..not.out1file) &
6077 write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
6078 " nfrag_back",nfrag_back
6080 read(inp,*) mset(iset)
6082 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
6084 if(me.eq.king.or..not.out1file) &
6085 write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
6086 ifrag(2,i,iset), qinfrag(i,iset)
6089 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
6091 if(me.eq.king.or..not.out1file) &
6092 write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
6093 ipair(2,i,iset), qinpair(i,iset)
6096 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
6097 wfrag_back(3,i,iset),&
6098 ifrag_back(1,i,iset),ifrag_back(2,i,iset)
6099 if(me.eq.king.or..not.out1file) &
6100 write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
6101 wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
6105 end subroutine read_fragments
6106 !-----------------------------------------------------------------------------
6108 !-----------------------------------------------------------------------------
6112 ! implicit real*8 (a-h,o-z)
6113 ! include 'DIMENSIONS'
6114 ! include 'COMMON.CSA'
6115 ! include 'COMMON.BANK'
6116 ! include 'COMMON.IOUNITS'
6118 ! integer :: ntf,ik,iw_pdb
6119 ! real(kind=8) :: var,ene
6121 open(icsa_in,file=csa_in,status="old",err=100)
6122 read(icsa_in,*) nconf
6123 read(icsa_in,*) jstart,jend
6124 read(icsa_in,*) nstmax
6125 read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
6126 read(icsa_in,*) nran0,nran1,irr
6127 read(icsa_in,*) nseed
6128 read(icsa_in,*) ntotal,cut1,cut2
6129 read(icsa_in,*) estop
6130 read(icsa_in,*) icmax,irestart
6131 read(icsa_in,*) ntbankm,dele,difcut
6132 read(icsa_in,*) iref,rmscut,pnccut
6133 read(icsa_in,*) ndiff
6140 end subroutine csa_read
6141 !-----------------------------------------------------------------------------
6142 subroutine initial_write
6145 ! implicit real*8 (a-h,o-z)
6146 ! include 'DIMENSIONS'
6147 ! include 'COMMON.CSA'
6148 ! include 'COMMON.BANK'
6149 ! include 'COMMON.IOUNITS'
6151 ! integer :: ntf,ik,iw_pdb
6152 ! real(kind=8) :: var,ene
6154 open(icsa_seed,file=csa_seed,status="unknown")
6155 write(icsa_seed,*) "seed"
6157 #if defined(AIX) || defined(PGI)
6158 open(icsa_history,file=csa_history,status="unknown",&
6161 open(icsa_history,file=csa_history,status="unknown",&
6164 write(icsa_history,*) nconf
6165 write(icsa_history,*) jstart,jend
6166 write(icsa_history,*) nstmax
6167 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
6168 write(icsa_history,*) nran0,nran1,irr
6169 write(icsa_history,*) nseed
6170 write(icsa_history,*) ntotal,cut1,cut2
6171 write(icsa_history,*) estop
6172 write(icsa_history,*) icmax,irestart
6173 write(icsa_history,*) ntbankm,dele,difcut
6174 write(icsa_history,*) iref,rmscut,pnccut
6175 write(icsa_history,*) ndiff
6177 write(icsa_history,*)
6180 open(icsa_bank1,file=csa_bank1,status="unknown")
6181 write(icsa_bank1,*) 0
6185 end subroutine initial_write
6186 !-----------------------------------------------------------------------------
6187 subroutine restart_write
6190 ! implicit real*8 (a-h,o-z)
6191 ! include 'DIMENSIONS'
6192 ! include 'COMMON.IOUNITS'
6193 ! include 'COMMON.CSA'
6194 ! include 'COMMON.BANK'
6196 ! integer :: ntf,ik,iw_pdb
6197 ! real(kind=8) :: var,ene
6199 #if defined(AIX) || defined(PGI)
6200 open(icsa_history,file=csa_history,position="append")
6202 open(icsa_history,file=csa_history,access="append")
6204 write(icsa_history,*)
6205 write(icsa_history,*) "This is restart"
6206 write(icsa_history,*)
6207 write(icsa_history,*) nconf
6208 write(icsa_history,*) jstart,jend
6209 write(icsa_history,*) nstmax
6210 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
6211 write(icsa_history,*) nran0,nran1,irr
6212 write(icsa_history,*) nseed
6213 write(icsa_history,*) ntotal,cut1,cut2
6214 write(icsa_history,*) estop
6215 write(icsa_history,*) icmax,irestart
6216 write(icsa_history,*) ntbankm,dele,difcut
6217 write(icsa_history,*) iref,rmscut,pnccut
6218 write(icsa_history,*) ndiff
6219 write(icsa_history,*)
6220 write(icsa_history,*) "irestart is: ", irestart
6222 write(icsa_history,*)
6226 end subroutine restart_write
6227 !-----------------------------------------------------------------------------
6229 !-----------------------------------------------------------------------------
6230 subroutine write_pdb(npdb,titelloc,ee)
6232 ! implicit real*8 (a-h,o-z)
6233 ! include 'DIMENSIONS'
6234 ! include 'COMMON.IOUNITS'
6235 character(len=50) :: titelloc1
6236 character*(*) :: titelloc
6237 character(len=3) :: zahl
6238 character(len=5) :: liczba5
6240 integer :: npdb !,ilen
6244 ! real(kind=8) :: var,ene
6248 if (npdb.lt.1000) then
6249 call numstr(npdb,zahl)
6250 open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
6252 if (npdb.lt.10000) then
6253 write(liczba5,'(i1,i4)') 0,npdb
6255 write(liczba5,'(i5)') npdb
6257 open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
6259 call pdbout(ee,titelloc1,ipdb)
6262 end subroutine write_pdb
6263 !-----------------------------------------------------------------------------
6265 !-----------------------------------------------------------------------------
6266 subroutine write_thread_summary
6267 ! Thread the sequence through a database of known structures
6268 use control_data, only: refstr
6270 use energy_data, only: n_ene_comp
6272 ! implicit real*8 (a-h,o-z)
6273 ! include 'DIMENSIONS'
6275 use MPI_data !include 'COMMON.INFO'
6278 ! include 'COMMON.CONTROL'
6279 ! include 'COMMON.CHAIN'
6280 ! include 'COMMON.DBASE'
6281 ! include 'COMMON.INTERACT'
6282 ! include 'COMMON.VAR'
6283 ! include 'COMMON.THREAD'
6284 ! include 'COMMON.FFIELD'
6285 ! include 'COMMON.SBRIDGE'
6286 ! include 'COMMON.HEADER'
6287 ! include 'COMMON.NAMES'
6288 ! include 'COMMON.IOUNITS'
6289 ! include 'COMMON.TIME1'
6291 integer,dimension(maxthread) :: ip
6292 real(kind=8),dimension(0:n_ene) :: energia
6294 integer :: i,j,ii,jj,ipj,ik,kk,ist
6295 real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
6297 write (iout,'(30x,a/)') &
6298 ' *********** Summary threading statistics ************'
6299 write (iout,'(a)') 'Initial energies:'
6300 write (iout,'(a4,2x,a12,14a14,3a8)') &
6301 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
6302 'RMSnat','NatCONT','NNCONT','RMS'
6303 ! Energy sort patterns
6308 enet=ener(n_ene-1,ip(i))
6311 if (ener(n_ene-1,ip(j)).lt.enet) then
6313 enet=ener(n_ene-1,ip(j))
6325 ist=nres_base(2,ii)+ipatt(2,i)
6327 energia(i)=ener0(kk,i)
6329 etot=ener0(n_ene_comp+1,i)
6330 rmsnat=ener0(n_ene_comp+2,i)
6331 rms=ener0(n_ene_comp+3,i)
6332 frac=ener0(n_ene_comp+4,i)
6333 frac_nn=ener0(n_ene_comp+5,i)
6336 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
6337 i,str_nam(ii),ist+1,&
6338 (energia(print_order(kk)),kk=1,nprint_ene),&
6339 etot,rmsnat,frac,frac_nn,rms
6341 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
6342 i,str_nam(ii),ist+1,&
6343 (energia(print_order(kk)),kk=1,nprint_ene),etot
6346 write (iout,'(//a)') 'Final energies:'
6347 write (iout,'(a4,2x,a12,17a14,3a8)') &
6348 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
6349 'RMSnat','NatCONT','NNCONT','RMS'
6353 ist=nres_base(2,ii)+ipatt(2,i)
6355 energia(kk)=ener(kk,ik)
6357 etot=ener(n_ene_comp+1,i)
6358 rmsnat=ener(n_ene_comp+2,i)
6359 rms=ener(n_ene_comp+3,i)
6360 frac=ener(n_ene_comp+4,i)
6361 frac_nn=ener(n_ene_comp+5,i)
6362 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
6363 i,str_nam(ii),ist+1,&
6364 (energia(print_order(kk)),kk=1,nprint_ene),&
6365 etot,rmsnat,frac,frac_nn,rms
6367 write (iout,'(/a/)') 'IEXAM array:'
6368 write (iout,'(i5)') nexcl
6370 write (iout,'(2i5)') iexam(1,i),iexam(2,i)
6372 write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
6373 'Max. time for threading step ',max_time_for_thread,&
6374 'Average time for threading step: ',ave_time_for_thread
6376 end subroutine write_thread_summary
6377 !-----------------------------------------------------------------------------
6378 subroutine write_stat_thread(ithread,ipattern,ist)
6380 use energy_data, only: n_ene_comp
6382 ! implicit real*8 (a-h,o-z)
6383 ! include "DIMENSIONS"
6384 ! include "COMMON.CONTROL"
6385 ! include "COMMON.IOUNITS"
6386 ! include "COMMON.THREAD"
6387 ! include "COMMON.FFIELD"
6388 ! include "COMMON.DBASE"
6389 ! include "COMMON.NAMES"
6390 real(kind=8),dimension(0:n_ene) :: energia
6392 integer :: ithread,ipattern,ist,i
6393 real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
6395 #if defined(AIX) || defined(PGI)
6396 open(istat,file=statname,position='append')
6398 open(istat,file=statname,access='append')
6401 energia(i)=ener(i,ithread)
6403 etot=ener(n_ene_comp+1,ithread)
6404 rmsnat=ener(n_ene_comp+2,ithread)
6405 rms=ener(n_ene_comp+3,ithread)
6406 frac=ener(n_ene_comp+4,ithread)
6407 frac_nn=ener(n_ene_comp+5,ithread)
6408 write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
6409 ithread,str_nam(ipattern),ist+1,&
6410 (energia(print_order(i)),i=1,nprint_ene),&
6411 etot,rmsnat,frac,frac_nn,rms
6414 end subroutine write_stat_thread
6415 !-----------------------------------------------------------------------------
6417 subroutine readpdb_template(k)
6418 ! Read the PDB file for read_constr_homology with read2sigma
6419 ! and convert the peptide geometry into virtual-chain geometry.
6421 ! include 'DIMENSIONS'
6422 ! include 'COMMON.LOCAL'
6423 ! include 'COMMON.VAR'
6424 ! include 'COMMON.CHAIN'
6425 ! include 'COMMON.INTERACT'
6426 ! include 'COMMON.IOUNITS'
6427 ! include 'COMMON.GEO'
6428 ! include 'COMMON.NAMES'
6429 ! include 'COMMON.CONTROL'
6430 ! include 'COMMON.FRAG'
6431 ! include 'COMMON.SETUP'
6432 use compare_data, only:nhfrag,nbfrag
6433 integer :: i,j,k,ibeg,ishift1,ires,iii,ires_old,ishift,ity, &
6435 logical lprn /.false./,fail
6436 real(kind=8), dimension (3):: e1,e2,e3
6437 real(kind=8) :: dcj,efree_temp
6441 real(kind=8), dimension (3,40) :: sccor
6443 integer, dimension (:), allocatable :: iterter
6444 if(.not.allocated(iterter))allocate(iterter(nres))
6451 write (2,*) "UNRES_PDB",unres_pdb
6459 read (ipdbin,'(a80)',end=10) card
6460 if (card(:3).eq.'END') then
6462 else if (card(:3).eq.'TER') then
6465 itype(ires_old-1,1)=ntyp1
6466 iterter(ires_old-1)=1
6467 #if defined(WHAM_RUN) || defined(CLUSTER)
6468 if (ires_old.lt.nres) then
6470 itype(ires_old,1)=ntyp1
6472 #if defined(WHAM_RUN) || defined(CLUSTER)
6476 ! write (iout,*) "Chain ended",ires,ishift,ires_old
6479 dc(j,ires)=sccor(j,iii)
6482 call sccenter(ires,iii,sccor)
6485 ! Fish out the ATOM cards.
6486 if (index(card(1:4),'ATOM').gt.0) then
6487 read (card(12:16),*) atom
6488 ! write (iout,*) "! ",atom," !",ires
6489 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
6490 read (card(23:26),*) ires
6491 read (card(18:20),'(a3)') res
6492 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
6493 ! & " ires_old",ires_old
6494 ! write (iout,*) "ishift",ishift," ishift1",ishift1
6495 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
6496 if (ires-ishift+ishift1.ne.ires_old) then
6497 ! Calculate the CM of the preceding residue.
6501 dc(j,ires_old)=sccor(j,iii)
6504 call sccenter(ires_old,iii,sccor)
6508 ! Start new residue.
6509 if (res.eq.'Cl-' .or. res.eq.'Na+') then
6512 else if (ibeg.eq.1) then
6513 ! write (iout,*) "BEG ires",ires
6515 if (res.ne.'GLY' .and. res.ne. 'ACE') then
6519 ires=ires-ishift+ishift1
6521 ! write (iout,*) "ishift",ishift," ires",ires,
6522 ! & " ires_old",ires_old
6523 ! write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
6525 else if (ibeg.eq.2) then
6527 ishift=-ires_old+ires-1
6529 ! write (iout,*) "New chain started",ires,ishift
6532 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
6533 ires=ires-ishift+ishift1
6536 if (res.eq.'ACE' .or. res.eq.'NHE') then
6539 itype(ires,1)=rescode(ires,res,0,1)
6542 ires=ires-ishift+ishift1
6544 ! write (iout,*) "ires_old",ires_old," ires",ires
6545 ! if (card(27:27).eq."A" .or. card(27:27).eq."B") then
6548 ! write (2,*) "ires",ires," res ",res," ity",ity
6549 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
6550 res.eq.'NHE'.and.atom(:2).eq.'HN') then
6551 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
6552 ! write (iout,*) "backbone ",atom ,ires,res, (c(j,ires),j=1,3)
6554 write (iout,'(2i3,2x,a,3f8.3)') &
6555 ires,itype(ires,1),res,(c(j,ires),j=1,3)
6559 sccor(j,iii)=c(j,ires)
6561 if (ishift.ne.0) then
6562 ires_ca=ires+ishift-ishift1
6566 ! write (*,*) card(23:27),ires,itype(ires)
6567 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and.&
6568 atom.ne.'N' .and. atom.ne.'C' .and.&
6569 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and.&
6570 atom.ne.'OXT' .and. atom(:2).ne.'3H') then
6571 ! write (iout,*) "sidechain ",atom
6573 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
6577 10 if(me.eq.king.or..not.out1file) &
6578 write (iout,'(a,i5)') ' Nres: ',ires
6579 ! Calculate dummy residue coordinates inside the "chain" of a multichain
6582 write(2,*) "tutaj",ires,nres
6584 ! write (iout,*) i,itype(i),itype(i+1)
6585 if (itype(i,1).eq.ntyp1.and.iterter(i).eq.1) then
6586 if (itype(i+1,1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
6587 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
6588 ! first is connected prevous chain (itype(i+1).eq.ntyp1)=true
6589 ! second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
6591 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
6592 call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
6599 c(j,i)=c(j,i-1)-1.9d0*e2(j)
6603 dcj=(c(j,i-2)-c(j,i-3))/2.0
6604 if (dcj.eq.0) dcj=1.23591524223
6609 else !itype(i+1).eq.ntyp1
6611 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
6612 call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
6619 c(j,i)=c(j,i+1)-1.9d0*e2(j)
6623 dcj=(c(j,i+3)-c(j,i+2))/2.0
6624 if (dcj.eq.0) dcj=1.23591524223
6629 endif !itype(i+1).eq.ntyp1
6630 endif !itype.eq.ntyp1
6632 ! Calculate the CM of the last side chain.
6635 dc(j,ires)=sccor(j,iii)
6638 call sccenter(ires,iii,sccor)
6642 if (itype(nres,1).ne.10) then
6646 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
6647 call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
6654 c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
6658 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
6659 if (dcj.eq.0) dcj=1.23591524223
6660 c(j,nres)=c(j,nres-1)+dcj
6661 c(j,2*nres)=c(j,nres)
6672 c(j,2*nres)=c(j,nres)
6674 if (itype(1,1).eq.ntyp1) then
6678 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
6679 call refsys(2,3,4,e1,e2,e3,fail)
6686 c(j,1)=c(j,2)-1.9d0*e2(j)
6690 dcj=(c(j,4)-c(j,3))/2.0
6696 ! Copy the coordinates to reference coordinates
6702 ! Calculate internal coordinates.
6703 if (out_template_coord) then
6704 write (iout,'(/a)') &
6705 "Cartesian coordinates of the reference structure"
6706 write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
6707 "Residue ","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
6709 write (iout,'(a3,1x,i4,3f8.3,5x,3f8.3)')&
6710 restyp(itype(ires,1),1),ires,(c(j,ires),j=1,3),&
6711 (c(j,ires+nres),j=1,3)
6714 ! Calculate internal coordinates.
6715 call int_from_cart(.true.,out_template_coord)
6716 call sc_loc_geom(.false.)
6718 thetaref(i)=theta(i)
6723 dc(j,i)=c(j,i+1)-c(j,i)
6724 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
6729 dc(j,i+nres)=c(j,i+nres)-c(j,i)
6730 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
6732 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
6733 ! & vbld_inv(i+nres)
6738 cref(j,i+nres,1)=c(j,i+nres)
6748 end subroutine readpdb_template
6749 !-----------------------------------------------------------------------------
6751 !-----------------------------------------------------------------------------
6752 end module io_config